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Abstract 

A  simple  dynamic  model  is  used  with  various  observations  to  provide  an  approx¬ 
imate  spectral  description  of  low  frequency  oceanic  variability.  Such  a  spectrum  has 
wide  application  in  oceanography,  including  the  optimal  design  of  observational  strat¬ 
egy  for  the  deployment  of  floats,  the  study  of  Lagrangian  statistics  and  the  estimate 
of  uncertainty  for  heat  content  and  mass  flux.  Analytic  formulas  for  the  frequency 
and  wavenumber  spectra  of  any  physical  variable,  and  for  the  cross  spectra  between 
any  two  different  variables  for  each  vertical  mode  of  the  simple  dynamic  model  are 
derived.  No  heat  transport  exists  in  the  model.  No  momentum  flux  exists  either  if 
the  energy  distribution  is  isotropic.  It  is  found  that  all  model  spectra  are  related  to 
each  other  through  the  frequency  and  wavenumber  spectrum  of  the  stream-function 
for  each  mode,  A),  where  (A;,  /)  represent  horizontal  wavenumbers,  w 

stands  for  frequency,  n  is  vertical  mode  number,  and  (<f>,  A)  are  latitude  and  longi¬ 
tude,  respectively.  Given  $(&,  l,u,  n,  <j>.  A),  any  model  spectrum  can  be  estimated.  In 
this  study,  an  inverse  problem  is  faced:  A)  is  unknown;  however,  some 

observational  spectra  are  available.  I  want  to  estimate  T ( k ,  l,u,  n ,  0,  A)  if  it  exists. 

Estimated  spectra  of  the  low  frequency  variability  are  derived  from  various  mea¬ 
surements.  (i)  The  vertical  structure  of  and  kinetic  energy  and  potential  energy  is 
inferred  from  current  meter  and  temperature  mooring  measurements,  respectively, 
(ii)  Satellite  altimetry  measurements  produce  the  geographic  distributions  of  surface 
kinetic  energy  magnitude  and  the  frequency  and  wavenumber  spectra  of  sea  surface 
height,  (iii)  XBT  measurements  yield  the  temperature  wavenumber  spectra  and  their 
depth  dependence,  (v)  Current  meter  and  temperature  mooring  measurements  pro¬ 
vide  the  frequency  spectra  of  horizontal  velocities  and  temperature. 

It  is  found  that  a  simple  form  for  A)  does  exist  and  an  analytical 

formula  for  a  geographically  varying  T ( A: ,  /,  to,  n,  (f>,  A)  is  constructed.  Only  the  energy 
magnitude  depends  on  location.  The  wavenumber  spectral  shape,  frequency  spec¬ 
tral  shape  and  vertical  mode  structure  are  universal.  This  study  shows  that  motion 
within  the  large-scale  low-frequency  spectral  band  is  primarily  governed  by  quasi- 
geostrophic  dynamics  and  all  observations  can  be  simplified  as  a  certain  function  of 
$(k,l,u,n,<f>,  A). 
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The  low  frequency  variability  is  a  broad-band  process  and  Rossby  waves  are  partic¬ 
ular  parts  of  it.  Although  they  are  an  incomplete  description  of  oceanic  variability  in 
the  North  Pacific,  real  oceanic  motions  with  energy  levels  varying  from  about  10-40% 
of  the  total  in  each  frequency  band  are  indistinguishable  from  the  simplest  theoretical 
Rossby  wave  description.  At  higher  latitudes,  as  the  linear  waves  slow,  they  disappear 
altogether.  Non-equatorial  latitudes  display  some  energy  with  frequencies  too  high 
for  consistency  with  linear  theory;  this  energy  produces  a  positive  bias  if  a  lumped 
average  westward  phase  speed  is  computed  for  all  the  motions  present. 

Thesis  Supervisor:  Carl  Wunsch 
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Chapter  1 


Introduction 


The  general  circulation  of  the  ocean  varies  on  all  time  and  space  scales.  One  of  the 
main  hindrances  to  understanding  the  climate  change  is  the  lack  of  a  quantitative 
description  of  the  natural  variability  of  the  ocean.  This  thesis  is  intended  to  produce 
a  quantitative  algebraic  spectral  description  of  low  frequency  oceanic  variability.  By 
the  expression  “low  frequency  variability” ,  I  mean  the  time-dependent  motions  with 
a  time  scale  longer  than  the  inertial  period  and  shorter  than  a  few  years,  and  a 
spatial  scale  ranging  from  tens  of  km  to  thousands  of  km.  Because  the  data  used  in 
this  study  typically  span  a  few  years,  “very  low  frequency  variability” — the  decadal 
and  even  longer-period  variability — is  not  considered  here.  Because  of  the  strong 
inhomogeneity  of  eddy  energy  level,  the  wavenumber  spectral  representation  is  not 
a  complete  one  for  the  basin  scale  in  the  ocean.  However,  in  the  open  ocean  away 
from  major  boundary  currents,  within  the  spatial  scale  of  about  1000  km,  the  energy 
level  does  not  change  much  with  location,  and  the  observed  wavenumber  spectra 
suggest  that  energy  of  time-dependent  motions  is  dominated  by  motions  with  spatial 
scales  of  a  few  hundreds  km.  Therefore,  the  wavenumber  spectral  representation  is 
still  very  useful  within  the  scale  of  about  1000  km.  Here  I  am  only  concerned  with 
mesoscale  eddies  in  the  regions  away  from  major  currents  and  the  equator.  Such  a 
spectral  description  has  broad  application.  It  allows  one  to  describe  the  oceanic  state 
concisely,  to  test  the  ocean  general  circulation  model  conveniently,  to  readily  compare 
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theory  with  observations  and  to  easily  design  observational  strategy. 

My  method  is  analogous  to  that  of  Garrett  and  Munk  (1972)  applied  to  internal 
waves.  Using  linear  dynamics  under  the  hypothesis  of  spatial  homogeneity,  horizontal 
isotropy  and  vertical  symmetry  of  the  wave  field,  Garrett  and  Munk  (1972)  patched 
together  a  universal  simple  algebraic  representation  of  the  distribution  of  internal 
wave  energy  in  wavenumber/frequency  space  in  the  deep  ocean,  which  has  become 
known  as  the  Garrett-Munk  spectrum.  Later  Garrett  and  Munk  (1975)  gave  an  im¬ 
proved  version.  The  task  for  low  frequency  variability  is  much  more  difficult  than  the 
case  for  internal  waves,  because  low  frequency  variability  exhibits  large  geographical 
variations  of  energy  level  and  strong  anisotropy. 

The  general  circulation  of  the  ocean  is  a  multi-dimensional  and  multi-variate  pro¬ 
cess.  In  the  last  few  decades,  thanks  to  the  advancement  of  measurement  technology, 
oceanographers  have  obtained  measurements  for  sea  surface  height,  temperature,  hor¬ 
izontal  velocities,  pressure,  etc.  at  different  places,  depths  and  times.  It  is  important 
to  find  out  whether  all  these  different  data  types  can  be  related  to  each  other.  If  so, 
this  large  amount  of  data  can  be  described  concisely.  The  fundamental  idea  of  my 
work  is  summarized  in  figure  1-1.  As  displayed  there,  both  the  data  and  model  provide 
estimates  of  the  true  oceanic  state  and  each  model  spectrum  is  related  to  the  oth¬ 
ers  through  the  frequency  and  wavenumber  spectrum  of  the  stream-function  for  each 
mode,  $(k,l,u,n,  0,  A),  where  ( k,l )  are  horizontal  wavenumbers;  u  is  frequency;  n 
represents  vertical  mode  number;  and  ( <j) ,  A)  stand  for  latitude  and  longitude,  respec¬ 
tively.  Given  $(k,  l,  u,  n,  <f>,  A),  each  model  spectrum  can  be  estimated  and  compared 
with  the  corresponding  spectrum  of  observations.  Here  an  inverse  problem  is  faced: 
<$(k,  /,  u>,  n ,  4>,  A)  is  unknown  and  some  observational  spectra  are  available;  I  want  to 
estimate  $(k,  l,  u>,  n,  (f>,  A)  from  the  observations  if  $(k,  l,  oj,  n,  (j),  A)  exists.  We  don’t 
know,  a  priori,  if  there  exists  a  simple  form  for  $>(k,l,u,n,(f),  A).  Because  neither 
the  data  nor  the  model  is  a  perfect  representation  of  the  true  time-dependent  mo¬ 
tions  in  the  ocean.  The  data  are  contaminated  by  instrument  noises  and  aliased  by 
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high  frequency  and  small  scale  motions  such  as  tides  and  internal  waves.  Moreover, 
the  available  data  are  short  time  series.  Thus,  the  spectrum  of  the  data  is  a  biased 
estimate  of  the  true  spectrum  of  the  underlying  oceanic  process.  There  are  many 
drawbacks  of  my  model  as  described  in  chapter  2.  In  my  simple  model,  it  is  assumed 
that  oceanic  processes  are  linear  and  homogeneous  in  space;  the  background  is  at 
rest;  there  is  no  forcing  and  the  depth  of  the  ocean  is  uniform.  None  of  these  assump¬ 
tions  is  true  for  the  actual  ocean.  Therefore,  it  is  unclear  whether  each  spectrum  of 
measurements  is  consistent  with  the  corresponding  model  spectrum.  If  they  are  con¬ 
sistent,  $(k,  l,  u,  n,  (j),  A)  will  exist  and  all  spectra  of  measurements  can  be  simplified 
as  a  certain  function  of  $(k,  l,  u,  n,  </>,  A). 

Assuming  that  the  energy  distribution  is  horizontally  isotropic,  Zang  (1998)  con¬ 
structed  a  simple  form  for  $(&,  /,  lj,  n,  <j),  A).  A  new  version  of  the  model  is  developed 
here  without  the  assumption  of  isotropy.  This  thesis  is  made  up  of  nine  chapters.  In 
chapter  2,  a  simple  linearized  model  for  low  frequency  variability  is  introduced  and 
the  solutions  of  the  model  are  derived.  Based  on  these  solutions,  analytic  formulas 
for  model  spectra  are  obtained  in  chapter  3.  In  chapter  4,  the  observed  properties 
of  low  frequency  oceanic  variability  are  presented.  The  emphasis  is  on  answering  the 
following  questions:  how  is  the  low  frequency  variability  energy  distributed  among 
horizontal  wavenumber,  direction,  frequency,  vertical  mode  and  geography?  and  does 
there  exist  some  kind  of  universal  spectral  description  for  low  frequency  variability  in 
the  ocean.  In  chapter  5,  a  simple  analytic  algebraic  representation  of  $(&,  l,  u,  n,  <f>,  A) 
is  constructed.  Application  of  the  spectrum  is  given  in  chapter  6.  In  chapter  7,  a 
directional  form  of  <$>(h,l,u,n,<j),  A)  is  produced,  which  can  differentiate  westward¬ 
going  energy  from  eastward-going  energy.  Chapter  8  presents  the  observed  dispersion 
relationship  for  North  Pacific  Rossby  wave  motions.  The  conclusion  and  discussion 
are  presented  in  chapter  9.  In  appendix  A,  the  relationship  between  Rossby  waves 
and  topography  is  studied. 
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Figure  1-1:  All  measurements  and  models  of  the  ocean  are  connected  to  each  other  through 
$(&,  l,u,  n,  0,  A).  Given  $(fc,  Z,a>,  n,  <j!>,  A),  one  can  calculate  all  model  spectra. 
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Chapter  2 


Dynamic  model  for  low  frequency 
motion 


2.1  The  governing  equations 

Away  from  the  equator  and  beneath  the  upper  mixed  layer,  the  time-dependent  mo¬ 
tions  of  the  continuously  stratified  ocean  can  be  described  by  the  following  standard 
linearized  equations  as  a  first-order  approximation  (e.g.,  Gill,  1982): 


du 

dt 

dv 

di 


dp 

dt 


1  dp 

-fv  = - 7p, 

Po  dx 

(2.1) 

,  ,  1  dp 

Po  dy 

(2.2) 

II 

1 

03 

*Nf  |^ 

1 

(2.3) 

+  wdp°{z)=0, 
dz 

(2.4) 

dv  dw 

+  di  +  al =  °- 

(2.5) 

Here  it  is  assumed  that  the  mean  velocities  are  zero,  i.e.,  u0  =  v0  =  w0  =  0.  Equations 
(2.1)  and  (2.2)  are  the  horizontal  momentum  equations;  equation  (2.3)  is  the  hydro¬ 
static  equation;  equation  (2.4)  is  the  density  conservation  equation  and  equation  (2.5) 
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is  the  continuity  equation.  The  variables  u,  v  and  w  are  the  perturbation  velocities;  p 
is  the  perturbation  density,  p  is  the  perturbation  pressure;  /  is  the  Coriolis  parameter 
and  po(z)  is  the  density  of  the  resting  ocean.  In  the  Boussinesq  approximation,  p0 
is  treated  as  constant  in  equations  (2.1)  and  (2.2).  The  dynamic  effect  of  salinity 
variations  has  been  simply  neglected  here.  Because  of  the  limitations  of  the  data,  the 
low  frequency  and  large  scale  variations  of  salinity  and  their  dynamic  effect  are  still 
unclear.  As  discussed  in  detail  by  Pedlosky  (1987),  equations  (2.1)  to  (2.5)  are  good 
representations  of  mesoscale  eddies  which  have  length  scale  of  hundreds  of  km.  For 
the  motions  with  horizontal  length  scale  comparable  to  or  greater  than  the  radius  of 
the  earth,  spherical  coordinates  should  be  used  instead  and  the  horizontal  variations 
of  the  basic  density  field  cannot  be  ignored. 

A  single  equation  in  p  can  be  obtained  from  equation  (2.1)  to  (2.5), 


d3p  d3p  d 2  2  d  1  d2p 

dx2dt  +  dy2dt  +  W  +  '  ’dz^dzdt 


+P[f 


2  3? 

dx 


where 


N2(z) 


g_dp 
Po  dz’ 


&P  2f &P 
dxdt 2  dydt 
(2.6) 


(2.7) 


and 


Separating  variables,  Eq.  (2.6)  is, 


(2.8) 


p(x,  y,  z. t) 


n~- f  oo 

Yj  Pn(x,y,Z,t ) 


n— 0 


7l— ~hoo 

Y  Pn(x,y,t)Fn(z), 


n— 0 


(2.9) 


a  sum  over  orthonormal  vertical  modes,  with  Fn  (z)  satisfying, 


d  (  1  dFn(z)\ 

dz  yA/^2(2:)  dz  ) 


(2.10) 
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where  r\  is  a  separation  constant.  The  horizontal  structure  is  governed  by, 


(A  +  f)l^L  +  **  _  +  A  +  stf2^  -  £2.  _  =  0 

''dt2  J  ^dx2dt  dy2dt  "Of  1  '  f«  1  J  0i  dxdt2  1  dydt’ 


(2.11) 


Sometimes  a  second  vertical  structure  function,  (z)  is  useful.  Define, 


G,n(-2:) 


1  dFn{z) 
N2(z )  dz 


(2.12) 


where, 

~  +  rnN2(z)Gn(z)  =  0.  (2.13) 

Eqs.  (2.10,  2.13)  are  readily  solved  numerically.  Analytic  solutions  to  the  vertical 
equation  are  available  for  a  few  forms  of  N(z),  including  i\T(z)=constant,  and  the 
exponential  profile, 

N(z)  =  N0ea\  (2.14) 

(Garrett  and  Munk,  1972),  which  more  closely  resembles  actual  buoyancy  profiles  in 
the  ocean.  Following  these  authors,  I  will  adopt  (2.14)  where  No  =  0.007  s_1  and 
a  =  0.001  m-1  in  dimensional  form,  as  a  zero-order  fit  below  about  1  km. 


2.2  Vertical  representation 


In  contrast  to  the  internal  wave  case,  observations  strongly  suggest  the  presence  of 
vertical  standing  modes  for  low  frequency  variability.  The  rigid-lid  upper  and  lower 
boundary  conditions  are  used  here: 


w(z)  =  0,  z  =  0,  —h. 


(2.15) 


Or, 


Gn(z)  =  0,  z  =  0,  —h. 


In  this  study,  I  assume  there  is  no  topography  and  take  h=4500  m. 


(2.16) 
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The  solution  to  equation  (2.13)  (Garrett  and  Munk,  1972)  is 


Gn{z)  =  AnJ0(—eaz)  +  BnY0(—eaz), 

(X  (X 


(2.17) 


where  J0(z ),  Y0  (z)  are  the  Bessel  function  of  the  first  kind.  Boundary  conditions 
require 

MU)Y0((o)  -  Jo(6>)roK-*)  =  0,  (2.18) 


and 


where, 


,  _  Jo(Norn/a) 

' n~Y0(N0rn/a )  : 


e  = 


N(,rn  , 

- e 

a 


&  = 


N0rn 

) 

a 


£-h  = 


r_ne-ah 

a 


(2.19) 


(2.20) 


Mode  No. 

Eigenvalue  ( rn ) 

gravity-wave  phase  speed  (c„) 

Equivalent  Depth  ( hn ) 

0 

0 

oo 

OO 

1 

0.402  (s/m) 

2.48  (m/s) 

63  (cm) 

2 

0.861  (s/m) 

1.16  (m/s) 

14  (cm) 

3 

1.319  (s/m) 

0.76  (m/s) 

6  (cm) 

Table  2-1:  Eigenvalue,  gravity- wave  phase  speed  and  equivalent  depth. 

The  first  four  eigenvalues,  rn,  are  listed  in  Table  2-1.  The  Rossby  radius  of  defor¬ 
mation  for  the  n-th  mode  can  be  obtained  from  rn  through 


Rn  =  (2.21) 

fr„ 

Chelton  et  al.  (1998)  computed  the  global  1°  x  1°  climatologies  of  the  first  baroclinic 
gravity  wave  phase  speed  c\  and  the  Rossby  radius  of  deformation  R\  from  climato¬ 
logical  average  temperature  and  salinity  profiles.  The  value  of  c\  ranges  from  about 
1.5  m/s  at  high  latitudes  to  3.0  m/s  near  the  equator.  The  geographical  variability 
of  R\  is  dominated  by  an  inverse  dependence  on  /.  They  also  discussed  the  effects 
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of  earth  rotation,  stratification,  and  water  depth  on  cn  and  Rn  in  terms  of  WKBJ 
approximation. 

Substitution  of  equation  (2.19)  into  (2.17)  yields 


Gn{z)  =  AnWn(z), 


(2.22) 


Wn(z)  =  J0(— eaz)  -  ^r!y0(^e“), 
a  Yo(N0rn/a)  a 


(2.23) 


1  dGn  _  AndWn(z ) 

rn\Z)  ~  o  —  o  7“  • 


dz  f2  dz 


Following  custom,  I  normalize  the  vertical  solution  by  setting 


(2.24) 


/V(z)dz  =  1. 

J  —h 


(2.25) 


Substitution  of  equation  (2.24)  into  equation  (2.25)  yields 


K  r°  ,dwn(z).2 


r £  7-/j  dz 


fdz  =  1. 


(2.26) 


By  using  the  properties  of  Bessel  function  (Abramowitz  and  Stegun,  1964),  I 


obtain 


r\dW„(z).2 


)2dz  =  abn 


(2.27) 


where 


bn  =  F  PUm  =  ^[y0"2(?o)  -  n-J(f-i)]-  (2.28) 

J<i-1  TV 


The  vertical  normalization  yields 


\/  abn 


(2.29) 
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The  normalized  vertical  eigenfunctions  are 


K(z)  = 


I  1  dWn(z) 

abv  dz 


-  N0rne° 


ac. 


1 


%(&)' 


JiK)],  (2.30) 


Gn{z)  = 


which  satisfy  the  orthogonal  conditions 


(2.31) 


[  Fm(z)F„(z)dz  =  Snm,  (2.32) 

J —h 

f  N2  (z)Gm(z)Gn(z)dz  =  r25nm.  (2.33) 

J  —h 

where  6nm  is  the  Kronecker  delta  function.  Note  that  the  units  of  Fn(z )  and  Gn(z )  are 
m-1/2  and  s2m~3/2,  respectively.  The  normalized  vertical  eigenfunctions  Fn(z )  and 
Gn(z )  are  plotted  in  figure  2-1  for  n=0  to  3.  An  important  property  in  figure  2-1  is  the 
near-surface  intensification  of  Fn(z)  for  n=l  to  3.  Because  the  vertical  structure  of  the 
horizontal  kinetic  energy  is  proportional  to  F2(z),  the  vertical  structure  of  baroclinic 
horizontal  kinetic  energy  will  correspondingly  show  a  near-surface  intensification. 


2.3  Horizontally  propagating  waves 

The  horizontal  structure  of  each  mode  is  represented  in  the  form  of  propagating  waves: 

/-f-OO  r+OO  r- l~oo 

/  /  p(k,  l,  w,  7i)el27dfc:E+  y-cA> dkdldu ,  (2.34) 

-oo  J  —  oo  J — oo 

where  (k,  l )  stand  for  horizontal  wavenumber;  and  uj  represents  frequency.  In  order 
to  compare  with  the  measurements  readily,  the  cyclical  frequencies  and  wavenumbers 
are  used  here.  The  units  of  cyclical  frequencies  and  wavenumbers  are  cycles  per  day 
and  cycles  per  kilometer,  respectively.  Substitution  of  equation  (2.34)  into  (2.11) 
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yields  the  dispersion  relation 


2mv(f2 — (2Ttu))2)[(2irk)2 +(2irl)2 +r2(f2 — (2ttlo)2)]+ (3([(f2 — (2iuj)2)2Trk+8n2  flu>]  =  0. 

(2.35) 

For  low  frequency,  2iroo  «  /,  equation  (2.35)  can  be  simplified  as 


cu[(27rA;)2  +  ( 2ttI )2  +  f2r2n }  +  J3k  =  0. 


(2.36) 


For  any  given  mode  at  a  fixed  latitude,  there  is  a  maximum  frequency  (shortest 
period)  for  planetary  waves  given  by 


t^n—max 


(3Rn 

47T 


(2.37) 


Recall  that  the  vertical  structure  of  the  pressure  for  each  mode  is  represented  by 
the  function  Fn(z).  Therefore,  the  full  solution  for  each  mode  in  pressure  is 

pn(x,  y,  z,  t)  =  rrrm, I,  u,  »)F„(z)e^<to^^>]clMWu>.  (2.38) 

J  — oo  J —o o  7—oo 


The  low  frequency  oceanic  variability  is  a  broad-band  process  and  consists  of 
Rossby  waves,  eddies,  etc.  Rossby  waves  are  particular  parts  of  it  in  that  the  frequency 
and  wavenumber  of  Rossby  waves  must  satisfy  the  dispersion  relation,  equation  (2.36). 
The  full  solution  for  Rossby  waves  for  each  mode  in  pressure  is 


r+c o  r- f-oo  r+oo  jjf £ 

pnr  (x,  y,  z,  t)  =  J  ^  J_JP(K  l ,  n)Fn{z)5{u-j +  j^) 


pk 


(2.39) 


where  5  is  the  delta  function. 

For  low  frequency  variability,  the  horizontal  momentum  equations  can  be  approx¬ 
imated  by  the  geostrophic  relations: 


fv  = 


1  dp 
Po  dx’ 


(2.40) 


V]dkdldu}, 
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J  “  -  o  ■ 

Po  dy 

So  the  complete  wave  solutions  for  the  horizontal  velocities  for  each  mode  are 


(2.41) 


/+OO  /'“fOC  r 

/  / 

-oo  */ — oo  - 


+°°r»2?rfp(A;>f>fa;>n)FB(g)l  n<kx+ly_, 


OO  P  — OO  —  OO 


Jt)dkdlduj,  (2.42) 


w(«,*, x,t)  =  (2.43) 

J—oo  J—oo  J  —  00  /9n  / 


The  solutions  for  the  density  and  vertical  velocity  can  be  derived  from  the  hydrostatic 
equation  and  the  density  conservation  equation.  The  complete  wave  solutions  for  the 
density  and  vertical  velocity  are 


Pn(x,y,z,t)  =  -  rrrim: 

J—oo  J— OO  J —00  I  (J 

(2.44) 

.(«, V, a. <)  =  /+“/+7+“  eW*Wdkdldw.  (2.45) 

V  — 00  J—oo  J—oo  I  On  1 


00  */  —00  */  —00 


p(M,a;>nOTz)GB(z)l  i2)r(fc»+i». 


One  frequently  used  variable,  the  vertical  displacement  £,  is  related  to  the  vertical 
velocity  through 

8  C 

m  =  m  ■  (2A6) 

Equation  (2.46)  yields  the  solution  for  each  mode  in  vertical  displacement 


_  f+°°  r+°°  r+00  f p(k,l,LO,n)Gn(z) 


/-f-oo  r+oo  r+00 

- 

-00  J—oo  J—OO 


ei2i;(,kx-Hy-wt)  dkdldu),  (2.47) 


The  vertical  displacement  is  not  directly  measured  in  the  real  ocean  but  is  inferred 
from  temperature  time  series  by 


C  {x,y,z,t) 


0{x,y,z,  t) 

d0o  /  dz 


(2.48) 


From  equation  (2.48),  the  solution  for  each  mode  in  temperature  can  be  derived 

/+7+7+“[ «.]  (2.49) 

J—oo  J—oo  J—oo  On  OZ 
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As  shown  in  equation  (2.48),  the  temperature  variations  in  my  model  are  attributed 
only  to  the  vertical  advection  of  the  mean  vertical  temperature  profile.  Therefore, 
my  model  cannot  represent  sea  surface  temperature  changes  because  of  the  rigid-lid 
approximation.  It  cannot  be  applied  to  the  mixed  layer  either.  In  the  mixed  layer, 
horizontal  advection,  mixing  and  atmospheric  thermal  forcing  play  an  important  role 
in  the  temperature  changes. 


2.4  Summary 

In  this  chapter,  a  simple  dynamic  model  is  introduced  and  analytic  solutions  to  the 
model  are  obtained.  For  simplification,  define 


•0( k ,  l,  to,  n ) 


Pof 


(2.50) 


which  represents  the  Fourier  transform  of  the  stream-function  for  each  mode.  If  the 
signal  is  real,  7p*(k,  l ,  u>,  n)  =  $(-k,  -l,  -w,  n),  where  *  stands  for  complex  conjugate. 
The  solutions  for  each  mode  can  be  summarized  as 


pn(x,  y,  z,  t)  =  rrr^K  i,  u>,  n)pa(k,  I,  u>,  z,  nY^^dkM w,  (2.51) 
7  —  00  7  —  00  7  —  00 


m 


wn(a:,y,z,i)  =  [  °°  f  °°  f  $(k,l,u,n)ua(k,l,u,  z,n)ei2*{kx+ly  Ut)dkdldu ,  (2.52) 
7  —  00  7  —00  7  —00 

(:r,  y,  z,  t)  =  [+°°  f+C°  [+°°'ip(k,  l,  w,  n)u«(fc,  l,  u>,  z,  n)e*’2»(**+I»-wt)  dkdlduj ,  (2.53) 

7—oo  7—oo  7—oo 

(*,»,*,«)  =  <2-54) 

7—oo  7—oo  7—oo 

*,(*,,,*,«)  =  nWn(y,^WM1V,.)<a(b4Mil*,  (2-55) 

7  —00  7  —  oo  7  —  00 

<„(*,»,*,«)  =  k™  f+~  [+C°i>(k,l,u>,n)  (2.56) 

7—oo  7—oo  7—oo 

»n(x,y,z,t)  =  [^[^r°j(k,l,U,n)9.{k,l,U,z,n)j (2.57) 
7  —  00  7  —  oo  7  —  co 
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where  I  have  defined 


pa(k,  l,  u,  z,  n)  =  pofFn(z),  (2.58) 

ua(k,l,u,  z,n)  =  —i2nlFn(z),  (2.59) 

va(k,l,co,  z,n)  =  i2nkFn(z),  (2.60) 

wa(k ,  l,  lo,  2,  n)  —  i2itu)fGn(z),  (2.61) 

pa(k,l,u,z,n)  =  ~^N2{z)Gn{z),  (2.62) 

9 

(a{k,  l,  U,  z,  n)  =  —fGn(z),  (2.63) 

df) 

0a(k,l,u;,z,n)  =  f-^Gn(z).  (2.64) 
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Chapter  3 


Spectra  of  the  model 


In  chapter  2,  a  simple  dynamic  model  is  introduced  and  the  solutions  to  it  are  derived. 
In  this  chapter,  I  will  obtain  the  analytic  formulas  for  the  frequency  and  wavenumber 
spectra  of  horizontal  velocities,  pressure,  sea  surface  height,  temperature,  etc.  for  each 
mode  based  on  the  solutions  obtained  in  chapter  2.  The  formula  for  the  frequency 
spectrum  of  the  range-averaged  temperature  from  tomographic  data  will  also  be 
derived.  The  frequency  and  horizontal  wavenumber  spectra  to  be  obtained  in  this 
chapter  are  for  each  mode.  In  chapter  5,  I  will  compare  the  model  spectrum  with 
the  corresponding  observations  by  adding  up  the  contributions  from  the  first  few 
important  modes. 

3.1  Covariance 

Let  h(x,y,  z,t,n)  represent  the  solution  to  the  model  for  the  n-th  mode,  where 
h(x,y,z,t,n)  is  any  one  of  the  dependent  variables  un(x,y,z,t),  vn(x,y,z,t),  etc., 
which  can  be  written  as 

/-foo  r-f~oo  /*+CO  „ 

/  /  rp(k,l,u,n)ha(k,l,co,  z,n)et2^kx+ly~u)t^dkdlduj,  (3.1) 

-co  J— oo  J ~oo 
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where  the  characteristic  factor  ha(kj,u,z,n)  correspondingly  stands  for  any  of  the 
variables  ua(k,  l, to,  z,  n ),  va(k ,  l, u,  z,  n),  etc. 

For  a  spatially  homogeneous  and  temporally  stationary  process, 

<i>(k,l,u,n)i>*(k',l',u,n)  >=  5(k'  —  k)8(l'  —  l)8(u}' -u)  <  $(k,  l,u,  n)\2  >  .  (3.2) 

where  the  angle  bracket  represents  ensemble  average.  Define 

${k,l,u,ri)  =<  \^(k,l,u,n)\2  >,  (3.3) 

which  is  the  three-dimensional  frequency  and  wavenumber  spectrum  of  the  stream- 
function  for  each  vertical  mode.  For  a  real  signal,  $(— k,  —  l,  —  u>, n)  =  $(k,l,u>,n). 

For  a  spatially  homogeneous  and  temporally  stationary  process,  the  auto-covariance 
function  of  h(x,  y,  z,  t,  n)  is 

/+oo  r+oo  r+ oo 

/  /  ha(k,l,u,z,n)h*a(k,l,u,z,n)$(k,l,u,n)e~l2w('kT:c+lry~u}r'>dkdldu 

-oo  J — oo  J — oo 

(3.4) 

Define  another  variable  s(x,  y,  z,  t,  n),  which  is  different  from  h(x,  y ,  z,  t,  n),  to  be 
the  solution  to  the  model,  then  s(x,y,z,  f,  n)  can  be  written  as 

/+oo  r+oo  r+oo  _ 

/  /  i>(k,l,LO,n)sa(k,l,uj,z,n)el2'K(-kx+ly~ult'>dkdldu j.  (3.5) 

-oo  J— oo  J  —  oo 

By  definition,  the  cross-covariance  function  between  h(x,  y ,  z,  i,  n)  and  s(x,  y,  z,  f,  n) 
for  a  homogeneous  and  stationary  process  can  be  written  as 


r+ oo  r+oo  r+oo 

Rhs(rx,ry,T,z,n)=  /  /  ha(k,l,co,  z,n)s*a(kj,u,  z,n)^(kJ,u,n)e~l27r<'krx+lry~lVT'>dkdldu. 

J — oo  J —oo  J— oo 
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3.2  Three-dimensional  spectra 

The  three-dimensional  spectra  for  u,  v,  etc.  for  each  mode  are 


^p,u,v,w,p,9,c{k,l,^,n,  z)  —  [\pa\2,\ua\2 ,\va\2,\wa\2,\pa\  j  |$a|  >  |Ca|  ]$(k,l,W,n).  (3.7) 


If  the  full  three-dimensional  spectrum  of  any  variable  for  each  mode  is  known,  one 
can  estimate  $(k,l,u,n)  through  equation  (3.7).  For  example,  given  the  three- 
dimensional  spectrum  of  pressure  for  each  mode,  can  be  obtained  through 


$(k,liLU,n) 


T p(k,l,w,n,  z) 
\pa(k,l,uj,n,z)\2' 


(3.8) 


The  reader  is  reminded  the  function  pa{k,  l,  u,  n,  z)  is  defined  in  equation  (2.58).  The 
full  three-dimensional  spectra  of  each  mode  are  very  hard  to  obtain  for  low  frequency 
oceanic  variability.  TOPEX/POSEIDON  altimeter  data  provide  nearly  simultaneous 
observations  of  global  sea  surface  height  every  9.91  days,  from  which  one  can  obtain 
the  three-dimensional  spectrum  of  sea  surface  elevation.  Pressure  perturbations  at 
sea  surface  are  related  to  sea  elevation  changes  through  (Wunsch  and  Stammer,  1998) 


P(x,  y,  z  =  0,  t)  =  gprj(x,  y,  t ) . 


(3.9) 


One  thus  can  infer  the  three-dimensional  spectrum  of  the  pressure  perturbations  at 
sea  surface.  However,  the  vertical  structure  of  pressure  is  unclear.  Besides  satellite 
measurements,  a  sparse  array  of  current  meter  and  acoustic  moorings  can  yield  esti¬ 
mates  of  directional  spectra.  Munk  (unpublished  manuscript,  1988)  showed  how  to 
get  a  handle  on  the  directional  wave  energy  distribution  by  interpreting  the  cross¬ 
spectra  of  travel  times  of  acoustic  transmissions. 


3.3  Two-dimensional  spectra 


The  k,  l-spectra 


The  two-dimensional  (k,  1)  spectra  for  u ,  v,  etc.  for  each  mode  are 

/-foo 

[l«a|2,  \va\2;  \Pa\2,  \&a\2 ,  |C a|2]$(M,  W,  n)dw.  (3.10) 

-oo 

The  k,  u-spectra 

The  two-dimensional  (k,  u>)  spectra  for  u,  v,  etc.  for  each  mode  are 

/+OO 

[K|2,  |v«| 2,  bol2,  \0a\2,  |Co| 2}$(k,l,u,n)dl.  (3.11) 

-oo 

The  l,  u)  -spectra 

The  two-dimensional  (1,  u>)  spectra  for  u,  v,  etc.  for  each  mode  are 

/  +  OO 

[K|2,  \va  |2,  \Pa\2,  \Qa\2,  \Va\2]$ {k ,  l ,  CO ,  Tl) dk .  (3.12) 

-oo 

Two-dimensional  spectra  in  the  ocean  can  be  obtained  from  a  few  repeated  XBT 
lines  and  alongtrack  altimeter  data. 


3.4  One-dimensional  spectra 


ui- spectra 

The  one-dimensional  frequency  spectra  for  u,  v,  etc.  for  each  mode  are 

/+OO  y+oo 

/  [W2,  kl2,  \Pa\2,  \0a\2,  \Ca\2Mk,l,u,n)dkdl.  (3.13) 

-OO  J —oo 

The  one-dimensional  frequency  spectra  of  u,  v  are  related  to  the  three-dimensional 
spectra  of  pressure  through 

T  „  -1  _  rp(k,l,ut,n,z)(2irk)2] 

- — - , -  -  dkdl  (3.14) 


k- spectra 
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The  one-dimensional  zonal-wavenumber  spectra  for  u,v,  etc.,  for  each  mode  are 


n+OO 

[kl2,  H2,  \Pa\2,  kl2,  \(a\2mk,l,u,n)dldu.  (3.15) 

-OO 


l-spectra 

The  one-dimensional  meridional-wavenumber  spectra  for  each  mode  are 


n-foo 

[K|2,  |ua|2,  k|2,  |(?a|2,  \(a\2]$(k,l,u),n)dkdu).  (3.16) 

-oo 


Given  the  three-dimensional  spectrum  of  pressure,  one  can  calculate  the  wavenum¬ 
ber  spectrum  of  velocity  through 

^  n  N  r+°°\rp(k,l,u>,n,z)(2nl)2  rp{kj,u>,n,z){2ivk)2] 

T .,&»,*)=// - JO - , - W -  dldu,  (3.17) 


T  /,  _.x_  f+°°[+°0f?P(k,l)v,n,z)(2irl)2  Tp(k,l,uJ,n,  z)(2nk)2~\ 

T .,(l,n,*)-/  / - to - , - 753 -  dkdu,.  (3.18) 


3.5  Horizontal  coherence 


For  two  points  at  the  same  depth  separated  horizontally  by  a  vector 


R  —  Xj  Xi  —  iXxt  ^y)i 


(3.19) 


the  cross  spectrum  between  h(x,y,  z,t,n)  and  s(x  +  rx,y  +  ry,z,t,n)  at  these  two 
points  is 


Ths{rx,ry,u,z,n)  =  f  f  ha(k,l,u),n,  z)s*a(k,l,u,n,  z)Q(k,  l,uj,n)e  *( kTx+lry) dkdl . 
J— oo  J— OO 

(3.20) 
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By  definition,  the  coherence  is 


/-too  r- |-oo 

/  has*a$(k,  l,oo,  n)e~i{-krx+lr^  dkdl 

^ hs\rx,ry,uj,z,n )  =  —j= . ~°°  .... 

I  r+oo  r+oo  r+oo  r+oo 

\  /  \ha\2${k,  l,co,n)dkdl  /  /  |sa|2$(£;,Z,cj,n)cZA;cZZ 

V  *'—00  7—00  7  —  00  7—00 

(3,21) 

If  the  horizontal  separation  is  zero,  the  cross  spectrum  becomes 


/-boo  r+oo 

/  ha(k,l,co,n,z)s*a(k,l,co,n,z)$(k,l,co,n)dkdl.  (3.22) 

-00  7—00 

Substitution  of  h  and  s  with  w,  v  and  9  into  the  above  equation  gives 


/+00  r+oo 

/  ua{k,l,uj,n,z)v*a{k,l,co,n,z)'^(k,l,co,n)dkdl,  (3.23) 

-00  7—oo 


/+00  r  +  oo 

/  ua(k,l,co,n,z)9l(k,l,to,n,z)$(k,l,co,n)dkdl,  (3.24) 

-00  7—oo 

/+00  r+oo 

/  va(k,l,co,n,z)9*a(k,l,u,n,z)<$>{k,l,co,n)dkdl.  (3.25) 

-OO  7  —  00 


Because  ua{k,l,u,n,z )  is  an  odd  function  of  Z;  va{k,l,co,n,  z)  is  an  odd  function 
of  k ;  0(&,  Z,  w,  n,  z)  is  independent  of  k;  and  Z  and  <&(£;,  Z,  a;,  n)  =  $(-£;,  -Z,  -or,  n),  it 
is  easy  to  prove  that 


++0, 0,  co,  z,  n)  =  Tvg( 0, 0,  co,  z,  n)  =  0, 


(3.26) 


cug{ 0, 0,  co,  z,  n )  =  cvg( 0, 0,  co,  z,  n)  =  0. 


(3.27) 


3.6  Energy  level 

kinetic  energy 

The  kinetic  energy  for  each  mode  at  depth  z  is 


Ek(n,z)  =  -[<  \un{x,y,z,t)\2  >  +  <  \vn(x,y,  z,t)\2  >] 
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=  \  [  [  [  {\ua\2  +  \va\2)$(k,l,LO,n)dkdldu),  (3.28) 

Z  j  —  GO  J —  OO  J  —  OO 

and  the  kinetic  energy  per  unit  surface  area  for  each  mode  is 

fO 

KE(n)  =  /  Ek(n,z)dz 

J —h 

1  r+oo  r+oo  r+oo  r  0 

=  -  /  /  /  /  (\ua\2  +  \va\2)$(k,l,u,n)dzdkdldu.  (3.29) 

Z  J— oo  J —oo  J —oo  J —h 

potential  energy 

The  potential  energy  per  unit  surface  area  for  each  mode  is 
PE(n)  =  \j_hN2iz)  <  1C n(x,y,z,t)\2  >  dz 

1  r+ oo  r+00  r+oo  rO 

=  -  /  /  /  iV2|Ca|2$(A:,  l,  lo,  n)dzdkdldu .  (3.30) 

2  J —oo  J —oo  J— oo  J—h 


temperature  variance 

Since  potential  energy  is  inferred  from  the  temperature  perturbation,  the  tempera¬ 
ture  variance  is  frequently  used  as  a  substitute  for  potential  energy.  The  temperature 
variance  for  each  mode  at  depth  z  is: 

/+oo  r+00  r+oo 

/  /  \Oa\2$(k,l,uj,n)dkdlduj.  (3.31) 

-oo  J —oo  J —oo 


total  energy 

The  total  energy  per  unit  surface  area  for  each  mode  is 


TE(n )  =  KE(n)  +  PE(n) 

1  r+oo  r+oo  r+oo  rO 

=  -  /  /  /  (|rta|2  +  \va\2  +  N2\Ca\2)$(k,  l,  u,n)dkdlduj. (3.32) 

Z  J —OO  J  —  OO  J  —  OO  J—h 
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3.7  Heat  and  momentum  transport 

By  definition,  the  heat  transport  of  each  mode  at  depth  2  can  be  written  as 


r-jr 00  r- f-00  00 

<  un6n  >=  Rue{0, 0, 0,  n,  z)  —  ua(k,  l,u>,  n,  z)9*a{k,  l,cv,  n,  z)$(k,  l,u,  n)dkdlduj 

J  —  OO  J  —  OO  J  —  OO 

/+OO 

T„e(0,0,o;,2,n)du;,  (3.33) 

-OO 


< 


r+oo  /•+00 

I’n^n  >=  Rvb{0,  0, 0,  n,  z)  =  /  /  /  va(k,l,u>,n,  z)6l(k,l,uj,n,  z)$(k,l,LU,n)dkdldu 

J —00  J —00  J —OO 

/+OO 

T„e(0, 0,  a;,  z,  n)duj.  (3.34) 

-OO 


Substitution  of  equation  (3.26)  into  equation  (3.33)  and  (3.34)  yields 


^  un9n  ^  v n0n  —  0. 


(3.35) 


There  is  no  heat  transport  for  the  simple  linear  dynamic  model  with  a  flat  bot¬ 
tom  because  two  components  of  horizontal  velocity  and  temperature  are  in  spatial 
quadrature. 

The  momentum  transport  is 


<  UnV-n  >  = 


c-t 00  r+oo  r  +  00 

>=  ic(0,0,0,n,z)  =  /  /  /  Ma(M)w)n,zK(M,w,n,z)$(M,w)n)Mdw 

J —OO  j—  OO  J —00 

f+ OO  r-fOO  r+OO 

=  -T2(z)  /  /  /  47r2M<h(A:,/,a;,n)Mda;.  (3.36) 

J—oo  j —OO  J —OO 


If  <!>(&,  Z,u;,n)  is  isotropic,  the  momentum  transport  <  unwn  >  is  zero  because  k  and 
l  are  odd  functions. 


3.8  Acoustic  tomography  data 

Most  of  the  data  types  we  use  measure  velocity  or  temperature,  etc.  at  a  point. 
Tomographic  data  differ  in  that  they  integrate  spatially  through  the  ocean  thus  fil- 
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tering  out  high  wavenumbers.  Assuming  that  a  preliminary  inversion  of  the  data  has 
been  done,  tomographic  observations  are  reduced  to  a  purely  horizontal  average  of 
temperature  at  fixed  depths  (Munk  et  ah,  1996). 

Suppose  the  acoustic  ray  path  is  along  the  east-west  direction  and  the  equivalent 
horizontal  tomographic  path  length  is  L.  Then  it  is  readily  shown  that  the  frequency 
spectral  density  of  the  tomographic  temperature  average  is 

/-boo  r+oo 

/  9a6*a$(k,l,u,n)W(k,L)dkdl.  (3.37) 

-oo  J—oo 


where 


W(k,L)  = 


sin2(nkL ) 
(nkL)2 


(3.38) 


represents  the  spatial  filtering. 

Equation  (3.37)  suggests  that  tomographic  measurements  filter  out  short  waves 
along  the  direction  of  acoustic  ray  path  and  have  no  effect  on  the  waves  propagating 
perpendicular  to  the  acoustic  ray  path.  Therefore,  one  can  estimate  the  directional 
property  of  the  energy  density  <£>(&;,  l,  co,  n )  by  comparing  the  frequency  spectra  of 
tomographic  data  with  identical  path  length  and  different  orientation. 


3.9  Summary 

Analytic  formulas  for  the  spectra  of  a  simple  model  are  derived  in  this  chapter.  The 
spectrum  for  each  mode  can  be  summarized  as 

(1)  The  three-dimensional  spectrum  for  p 

T p{k,  l, a;,  n,  z)  =  plf2F%(z)$(k,  l,  w,  n),  (3.39) 

(2)  The  three-dimensional  spectrum  for  rj 

Tv(k,  l, o;,  n)  =  ^Fn(z  =  0 )$(&,  l, o;,  n),  (3.40) 
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(3)  The  three-dimensional  spectrum  for  u  and  v 


TUjV(k,l,uj,n,z)  =  [4-K2l2,4iT2k2]FZ(z)$(k,l,u,n), 


(3.41) 


(4)  The  three-dimensional  spectrum  for  9  and  £ 


T otC(k,l,u),n,z)=  f2G2n{z){—)2,f2G2n{z)  $(k,l,uj,n), 


(5)  The  three-dimensional  spectrum  for  w  and  p 


(6)  The  zonal-wavenumber/frequency  spectrum  for  p 


(3.42) 


T  W!P(k,l,u,n,z)  =  [47 i2uj2f2G2n(z),—  —  §(k,l,u,n),  (3.43) 


/■+ o° 

T p(k,uj,n,z)  =  p20f2F2(z)  $(£;,  l,  uj,  n)dl,  (3.44) 

J  —  OO 


(7)  The  meridional- wavenumber/frequency  spectrum  for  p 


r-\-oo 

Tp(l,  uj,  n,  z)  =  plf2F2{z )  /  $(k,  l ,  u,  n)dk, 

J  —  OO 


(3.45) 


(8)  The  frequency  spectrum  for  u ,  v,  p  and  9 


n-fOO 

4 tt212F2(z),  4tt 2k2F2(z),p2f2F2(z ),  f2G2n(z){—±)2  l,u ,  n)dkdl 

-00  nz 


(9)  The  frequency  spectrum  for  horizontally  averaged  temperature 


0+0°  >)( y 

f  G2n{z)(^-)2w{k,  L)${k,  l,u,  n)dkdl, 

-00  oz 


(10)  The  zonal- wavenumber  spectrum  for  u ,  v,  p  and  0 


(3.46) 


(3.47) 


,,r,.(k,n,z)  =  4»Vf ’(2),  ^(z),  /^(z)  A!1  i(k,l,u,n)dldu 


(11)  The  meridional-wavenumber  spectrum  for  u ,  v,  p  and  6 


TU}V>p<g(l,n,z)  =  YK2l2Fl{z),An2k2Fl(z),plf2Fl{z),f2G2n{z){-^Y  $(k,l,w,n)dkdu>, 

(3.49) 

(12)  The  kinetic  energy  for  each  mode  at  depth  2 


/*+ oo  f+ oo  r+oo 

Ek(n,  z)  =  F2(z)  /  /  /  2n2(k2  +  l2)$(k,  l,u,n)dkdlduj,  (3.50) 

2— oo  2  — OO  J— oo 


(13)  The  kinetic  energy  per  unit  surface  area  for  each  mode 


/-foo  r+oo  r+ oo 

/  /  2-K2(k2  +  l2)$(k,l,u,n)dkdlduj,  (3.51) 

-oo  2 — oo  2 — oo 


(14)  The  potential  energy  per  unit  surface  area  for  each  mode 


/+CO  r  +  OO  r+ OO  f  2^.2 

/  /  :^-iL$(A:,  l,  u,  n)dkdlduj, 

-oo  2 — oo  2— oo  2 


(3.52) 


(15)  The  total  energy  per  unit  surface  area  for  each  mode 


/-f OO  r+CO  r+OO 

/  /  [2tt2(*;2  +  Z2)  + +- ^(fc,l,w,n)dMdu,  (3.53) 

-oo  2  —  00  2  —  00  2 


(16)  The  ratio  of  kinetic  energy  to  potential  energy  per  unit  surface  area  for  each 


mode 


/-TOO  rH-OO  r+OO 

J  J  2'n2(k2  +  l2)$(k,  l,  a),  n)dkdldu) 
PE{n)  j+°°  r+°°  [+00(fr2j2)<t>(kJ,LU,n)dkdldu; 

2—oo  2~oo  2—oo 


(17)  The  heat  flux  for  each  mode 


(3.54) 


<  ^7+71  ^  ^ —  Oj 


(3.55) 


(18)  The  momentum  flux  for  each  mode 


/+oo  r+oo  r+oo 

/  /  47 r2kl$(k,l,(jj,n)dkdlduj.  (3.56) 

-OO  2—  00  2—00 


G  11‘n  Vji  ~ 
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If  $(k,l  ,  u>,  n )  is  isotropic,  the  momentum  transport  <  unvn  >  is  zero. 

The  general  circulation  of  the  ocean  is  a  multi-dimensional  and  multi-variate  pro¬ 
cess.  In  spite  of  its  complexity,  the  above  equations  suggest  that  different  spectra 
for  different  variables  are  connected  to  each  other  through  the  intermediate  variable 
which  is  the  three-dimensional  spectrum  of  the  stream-function  for  each 
mode.  Once  <b(k,l,cj,n)  is  known,  one  can  estimate  the  frequency  and  wavenumber 
spectra  for  any  variable  at  any  depth.  Accordingly,  the  correlation  functions,  time 
and  space  scales  of  temperature,  velocities,  etc.  can  be  obtained. 
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Chapter  4 


Observed  spectra  of  low  frequency 
oceanic  variability 


The  observed  spectra  of  low  frequency  variability  in  the  ocean  will  be  reviewed  in  this 
chapter.  “Low  frequency  variability”  is  defined  as  the  time-dependent  motions  with 
time  scales  longer  than  the  inertial  period  and  shorter  than  a  few  years  and  spatial 
scales  between  about  10  km  and  1000  km.  Emphases  are  put  on  answering  the  fol¬ 
lowing  questions:  How  does  the  variability  of  sea  surface  height,  horizontal  velocities 
and  temperature  depend  on  horizontal  wavenumbers,  frequency,  direction,  vertical 
mode  and  geography?  Which  characteristic  features  of  low  frequency  variability  are 
independent  of  geography?  Which  ones  are  geography-dependent,  and  how  do  they 
depend  on  geography?  I  will  summarize  the  various  observations  and  give  a  zero 
order  description  of  the  low  frequency  oceanic  variability. 

4.1  Horizontal  inhomogeneity 

The  most  striking  property  of  low  frequency  oceanic  variability  is  that  kinetic  energy 
varies  strongly  with  position.  This  strong  spatial  variation  implies  that  a  spectral 
representation  can  not  be  a  complete  one  because  it  vitiates  the  necessary  condition 
for  the  simple  spectral  representation:  spatial  homogeneity.  However,  if  the  kinetic 
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energy  does  not  change  much  within  the  scale  where  the  dominant  energy  of  the 
spectral  representation  occurs,  the  spectrum  is  still  a  useful  representation. 

The  spatial  variation  inferred  from  altimetry  and  current  meter  moorings  is  used 
as  first  approximation.  Figure  4-1. a  (taken  from  Wunsch,  1997)  shows  surface  kinetic 
energy  from  the  current  meters  in  the  North  Pacific;  figure  4-1. b  (taken  from  Stammer, 
1997a)  shows  the  surface  kinetic  energy  from  the  altimeter  in  the  North  Pacific;  and 
figure  4-l.c  shows  the  surface  kinetic  energy  from  the  empirical  formula,  equation 
(4.1).  Figure  4-2  is  the  same  as  figure  4-1  except  in  the  North  Atlantic.  There  is 
rough  agreement  between  the  two  estimates. 

The  zonal  averages  between  0°  and  360°E  of  eddy  kinetic  energy  KE  and  sea 
surface  slope  Ksi  =  KEsin2(f)  are  provided  as  a  function  of  latitude  in  figure  4-3  (taken 
from  Stammer  1997a).  The  zonal  averaged  KE  decreases  from  maximum  amplitude 
near  the  equator  to  minimum  amplitude  in  high  latitudes.  In  terms  of  Ksi,  values 
remain  almost  constant  in  the  low  latitudes  between  25°S  and  25°N.  This  implies  that 
equatorwards  of  25°  the  zonally  averaged  KE  is  approximately  in  inverse  proportion 
to  sin20. 

According  to  the  figures  4-1. a  and  4-1. b,  the  surface  eddy  kinetic  energy  in  the 
North  Pacific  can  be  regarded  as  being  composed  of  four  parts:  (1)  the  background 
part,  (2)  the  low  latitude  part  (south  of  the  energetic  currents)  where  KE  oc  (l/sin20), 
(3)  the  high  energy  source  with  a  center  around  (35°N,  150°E),  and  (4)  the  low  energy 
area  in  the  north  North  Pacific.  The  surface  eddy  kinetic  energy  in  the  North  Pacific 
is  expressed  as 


■Efcp(0,  A,  z  —  0)  —  30  + 


32 

sin20 


+  1000exp{  — [ 


(A-150)2  (0  —  35) 


900 


+ 


50 


-]} 


— 80exp{  — [ 


(A-190)2  ,  (0-42) 


1600 


+ 


200 


-]}, 


(4.1) 


where  A  is  the  longitude  from  0°  to  360°E,  and  0  is  the  latitude  from  10°  to  60°N. 
The  units  of  Ekp  are  cm2/s2. 
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log,0(KE)  (cm2/s!)  from  current  meters 


120°  150°  180°  210°  240°  270  °E 


T/P  log)0(KE)  (cm!/s!)  from  repeat  8-117 

60°N 


30° 


0.5  1.0  1.5  2.0  2.5  3.0  3.5 


Figure  4-1:  Surface  eddy  kinetic  energy  in  the  region  of  North  Pacific,  (a)  from 
the  current  meters  (Wunsch,  1997)  (b)  from  TOPEX/POSEIDON  data  (Stammer, 
1997a)  (c)  from  the  empirical  formula. 
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The  surface  eddy  kinetic  energy  in  the  North  Atlantic  is  expressed  as 


Eka(<l),\,z  =  0)  =  50  + 


sin  z(f> 


400 


80 


»>«P{-[^!  +  M],  -  160exp{-(^#>£  +  WiA]},  (4.2) 


2000 


200 


900 


50 


in  units  of  cm2/s2.  The  results  from  equations  (4.1)  and  (4.2)  are  drawn  in  figures 
4-1. c  and  4-2.c,  respectively.  As  shown  in  figures  4-1  and  4-2,  equations  (4.1)  and 
(4.2)  are  reasonable  fits  to  the  general  pattern  of  the  corresponding  observations. 

4.2  A  first-guess  k ,  /,  to  spectral  form 

The  spatial  inhomogeneity  of  oceanic  kinetic  energy  shows  that  there  is  no  universal 
form  $(k,l,u,n)  for  the  low  frequency  variability.  We  have  to  modify  the  univer¬ 
sal  form  $(k,l,u,ri)  to  a  regional  form  u),  n,  (j>,  A)  so  that  the  model  spec¬ 

trum  can  fit  the  corresponding  observations  at  different  locations.  Suppose  that 
A)  exists  and  a  convenient  representation  of  $(k,l,u),n,<j>,\)  for  each 

mode  is 

$(fc,  l, w,  n,  <f>,  A)  =  Bn(k)Cn{l)Dn(Lu)E0(n)I((f>,  A).  (4.3) 

Here  Bn(k )  and  Cn(l )  stand  for  the  zonal-  and  meridional- wavenumber  spectral  shape, 
respectively;  Dn(co)  represents  the  frequency  spectral  shape;  E0(n)  is  a  constant  asso¬ 
ciated  with  each  mode  which  will  determine  how  the  energy  is  divided  among  vertical 
modes;  and  I ((f),  A)  is  the  spatial  function  that  represents  how  the  energy  level  de¬ 
pends  on  space.  As  shown  in  figure  4-1  and  4-2,  in  the  areas  away  from  major  currents 
and  within  the  horizontal  scale  of  about  1000  km,  the  energy  level  can  be  roughly 
treated  as  constant,  and  thus  wavenumber  spectrum  is  a  useful  representation.  In 
equation  (4.3),  it  is  assumed  that  frequency  spectral  shape,  wavenumber  spectral 
shape,  vertical  structure  and  energy  level  are  separable,  which  is  the  simplest  case  for 
$(&,  l, u,  n,  (f),  A).  However,  as  will  be  seen  later,  the  separated  form  in  equation  (4.3) 


40 


cannot  distinguish  westward-going  energy  from  eastward-going  one.  A  sophisticated 
frequency  and  wavenumber  coupled  form  is  proposed  in  chapter  7  to  solve  this  prob¬ 
lem.  In  the  next  section,  I  will  decide  on  the  form  of  Bn(k),  Cn(l),  Dn(u),  and  E0(n ) 
according  to  different  observations. 

4.3  Vertical  structure  of  kinetic  energy  and  poten¬ 
tial  energy 

Schmitz  (1978,  1988)  studied  the  vertical  distribution  of  kinetic  energy  at  different 
sites  in  the  Atlantic  and  Pacific.  He  concluded  that  the  vertical  profile  for  kinetic 
energy  is  independent  of  geography  as  a  first  approximation  across  the  entire  mid¬ 
latitude  band.  Eddy  kinetic  energy  Ke  dropped  exponentially  from  the  surface  to  a 
depth  of  about  1.2  km,  then  remained  almost  constant  within  the  abyss. 

Wunsch  (1997)  systematically  studied  the  partition  of  the  kinetic  energy  of  time- 
varying  motions  amongst  the  dynamical  modes  throughout  the  water  column.  He 
concluded  that  in  the  open  North  Pacific,  the  water  column  average  kinetic  energy 
was  about  35%  in  the  barotropic  mode  and  about  55%  in  the  first  baroclinic  mode. 
The  North  Atlantic  contained  on  average  about  40%  in  the  barotropic  and  50%  in 
the  first  baroclinic  mode  in  the  middle  of  the  ocean.  Strong  deviations  from  these 
values  occur  near  the  Gulf  Stream  and  near  the  equator,  and  no  attempt  to  model 
them  in  those  regions  will  be  made  here.  The  results  by  Wunsch  and  Schmitz  are 
consistent.  Wunsch’s  result  shows  that  in  the  open  ocean  horizontal  kinetic  energy 
is  dominated  by  the  barotropic  and  first  baroclinic  mode  and  that  the  partitioning 
among  barotropic  and  first  baroclinic  mode  is  roughly  equal.  Because  the  first  baro¬ 
clinic  mode  is  surface  intensified,  horizontal  kinetic  energy  will  be  dominated  by  the 
first  baroclinic  mode  in  the  upper  ocean;  thus  it  will  decay  nearly  exponentially  in  the 
main  thermocline.  Within  the  deep  water,  horizontal  kinetic  energy  remains  almost 
constant  because  the  deep  ocean  is  dominated  by  the  barotropic  mode. 

Wunsch  (1999b)  inferred  the  vertical  displacement  of  an  isopycnal  from  tempera- 
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ture  measurements  and  investigated  the  vertical  structure  of  potential  energy  in  the 
North  Atlantic.  He  found  that  the  water  column  average  potential  energy  was  about 
30-40%  in  the  first  baroclinic  mode.  The  ratio  of  kinetic  energy  in  the  first  baroclinic 
mode  to  potential  energy  was  about  0.2-0. 4  in  the  mid-latitude  away  from  the  west¬ 
ern  boundary.  Reader  beware  that  the  results  for  the  vertical  structure  of  potential 
energy  are  semi-quantitative,  because  the  results  are  sensitive  to  the  values  of  39/ dz 
whose  reliability  is  not  clear  at  all,  and  the  number  of  useful  records  on  each  mooring 
is  very  limited. 

4.4  Observed  three-dimensional  spectra  of  SSH 

TOPEX/POSEIDON  satellite  altimeter  data  provide  nearly  simultaneous  measure¬ 
ments  of  global  sea  surface  height  every  10  days,  which  provides  estimates  of  the 
frequency /wavenumber  spectra  of  sea  surface  height.  The  spectra  were  first  derived 
by  Wunsch  and  Stammer  (1995)  from  the  spherical  harmonic  analysis  of  the  2°  x  2° 
gridded  data.  Subsequently,  Stammer  (1997a)  obtained  the  spectra  from  the  Fourier 
analysis  of  alongtrack  data.  The  spectra  obtained  from  these  two  different  methods 
are  quite  similar  (Wunsch  and  Stammer,  1995).  In  order  to  readily  compare  with 
model  results,  I  will  derive  the  spectra  from  the  Fourier  analysis  of  gridded  data 
to  see  how  the  spectra  depend  on  zonal-wavenumber,  meridional-wavenumber  and 
frequency.  Moreover,  the  spectra  produced  by  Wunsch  and  Stammer  (1995)  and 
Stammer  (1997a)  are  scalar  in  that  neither  Wunsch  and  Stammer  (1995)  nor  Stam¬ 
mer  (1997a)  studied  the  directional  property  of  the  spectrum.  In  this  chapter,  I  will 
provide  an  estimate  of  what  percentage  energy  in  the  sea  surface  height  propagates 
to  the  east  and  to  the  west.  The  value  used  here  is  gridded  with  a  spatial  resolution 
of  0.25°  x  0.25°  and  a  temporal  resolution  of  10  days.  These  numbers  were  produced 
by  Le  Traon  et  al.  (1998)  through  combining  TOPEX/POSEIDON  and  ERS-1/2 
measurements.  Because  the  low  frequency  oceanic  variability  strongly  depends  on 
geography  as  described  in  section  4.1,  I  will  study  the  three-dimensional  spectrum 
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region  by  region.  As  a  representative  result,  I  analyze  the  three-dimensional  spectra 
of  sea  surface  height  anomaly  in  the  area  of  25°  —  50°N  and  195°  —  225°E  in  North 
Pacific.  The  size  of  the  area  is  chosen  as  a  result  of  a  compromise  between  resolution 
and  homogeneity.  Because  the  meridional  range  of  the  area  is  so  small,  compared 
with  the  radius  of  the  earth,  this  area  can  be  assumed  to  be  in  a  horizontal  plane; 
that  is,  the  distance  of  one  degree  along  the  zonal  direction  is  independent  of  latitude. 
Defining  the  space-time  Fourier  transform  of  sea  surface  height  anomaly  r](x,  y,  t)  as 
fj(k,l,uj),  then 

/+OO  /»+ OO  Z'  +  OO 

/  /  r](x,y,t)el2^kx+ly-^dxdydt.  (4.4) 

-oo  J— oo  J — OO 

By  definition,  the  three-dimensional  spectrum  can  be  written  as 

~£r)(k,l,(jj)  — <  \fj(k,l,u)\2  >  .  (4.5) 

Because  the  data  are  real,  the  three-dimensional  spectrum  is  symmetrical  in  the  sense 
that  T v(k,  l,  oj)  =  Tjj( — k,  —l,  —u>).  Hereafter,  the  discussion  is  limited  to  positive  fre¬ 
quencies  only,  though  k  and  l  have  both  signs.  To  obtain  some  statistical  stability, 
the  spectrum  is  averaged  over  three  neighboring  frequencies.  This  produces  an  ap¬ 
proximate  6  degrees  of  freedom.  Figure  4-4  shows  the  smoothed  three-dimensional 
spectrum  at  periods  longer  than  60  days  plotted  as  log[Tv(k,  0/T7?-max(^)  01>  where 
Y,j_max(&,  0  is  the  peak  of  the  three-dimensional  spectrum  in  each  frequency  band. 
The  most  conspicuous  point  of  figure  4-4  is  that,  in  each  frequency  band,  energy  is 
dominated  by  long  wavelengths.  For  periods  longer  than  90  days,  the  westward-going 
energy  is  greater  than  eastward-going  energy  over  long  wavelengths.  In  each  frequency 
band,  the  energy  distribution  tends  to  be  isotropic  at  wavelengths  shorter  than  200 
km.  It  is  unclear  whether  this  isotropic  property  is  real  or  if  it  is  the  artificial  result  of 
the  objective  mapping.  The  spatial  correlated  function  used  to  produce  the  gridded 
data  is  isotropic  and  has  a  correlation  scale  about  200  km  in  the  area  we  studied  (Le 
Traon,  et  ah,  1998). 
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Figure  4-4:  The  smoothed  three-dimensional  spectra  at  periods  longer  than  59 
days.  The  spectra  at  shorter  periods  are  similar  except  that  the  energy  distribution 
is  nearly  isotropic. 
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4.5  Observed  two-dimensional  spectra  of  SSH 


As  seen  above,  the  three-dimensional  spectrum  is  too  complex  to  describe  in  an 
analytical  form.  It  can  be  simplified  as  two-dimensional  by  integrating  over  the 
meridional- wavenumber  l  or  the  zonal- wavenumber  k. 
k  —  u  spectrum 

By  definition,  the  zonal- wavenumber/frequency  spectrum  is 

/+oo 

Tv(k,l,u)dl.  (4.6) 

-OO 

Tv(k,u)  can  be  broken  down  into  three  parts:  the  eastward-going  T v{k  >  0,u  >  0), 
westward-going  T v(k  <  0,to  >  0),  and  standing  part  Tr](k  =  0,u  >  0)  (figure  4- 
5).  The  total  energy  over  the  positive  frequencies  is  14.2  cm2,  among  which  1.8 
cm2  is  contained  in  the  eastward-going,  4.1  cm2  in  the  westward-going  and  8.4  cm2 
in  the  standing  part.  Roughly  speaking,  the  energy  of  eastward-going,  westward¬ 
going,  and  standing  motions  is  in  the  ratio  of  1:2:4.  As  expected,  the  energy  of 
westward-going  motions  is  greater  than  that  of  eastward-going  motions.  The  spectra 
of  the  eastward-going  and  westward-going  motions  are  red  in  both  frequency  and 
wavenumber.  A  peak  occurs  near  the  period  of  60  days,  which  is  due  to  the  aliasing 
of  the  tidal  residuals.  Although  the  total  energy  of  the  westward-going  part  is  more 
than  twice  that  of  eastward-going  part,  it  is  hard  to  tell  this  difference  from  figure  4-5. 
To  further  elucidate  this  difference,  the  two-dimensional  spectrum  is  integrated  over 
zonal-wavenumber  or  frequency  to  generate  the  one-dimensional  spectrum,  which  will 
be  discussed  later.  More  than  half  of  the  total  energy  is  contained  in  the  standing 
part.  The  standing  part,  which  is  red,  is  a  function  only  of  frequency.  The  most 
conspicuous  peak  in  figure  4-5  is  associated  with  the  annual  cycle,  and  the  ridge  near 
the  period  of  60  days  again  results  from  the  aliasing  of  the  tides. 

I  —  uj  spectrum 

Similarly,  the  two-dimensional  meridional- wavenumber/frequency  spectrum  T v(l,  u) 
can  be  obtained  by  integrating  the  three-dimensional  spectrum,  Tv(k,l,ou),  over  the 
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zonal- wavenumber.  The  spectra  of  the  northward-going  Tv(l  >  0,  u)  >  0),  southward¬ 
going  T v(l  <  0,6;  >  0)  and  standing  part  T v(l  =  0,  co  >  0)  are  shown  in  figure 
4-6.  The  total  energy  over  the  positive  frequencies  is  14.2  cm2,  among  which  2.9 
cm2  is  contained  in  the  northward-going,  3.9  cm2  in  the  southward-going  and  7.4  cm2 
in  the  standing  motions.  The  total  energy  of  northward-going  and  southward-going 
motions  is  approximately  the  same  as  the  energy  of  standing  motions.  Although 
there  is  more  energy  going-southward,  the  difference  between  northward- going  and 
southward-going  energy  is  much  smaller  than  the  difference  between  the  eastward¬ 
going  and  westward- going  energy. 

4.6  Observed  one-dimensional  spectra  of  SSH 

As  discussed  in  last  section,  two-dimensional  spectra  are  still  fairly  complex.  They 
can  be  further  simplified  into  one-dimensional  forms. 

The  zonal-wavenumber  spectra 

The  integration  of  zonal-wavenumber/frequency  spectrum  of  the  eastward-going 
motions  over  frequency  produces  the  zonal-wavenumber  spectrum  of  eastward-going 
motions, 

r+oo 

T r,e(k)=JQ  Tv(k  >  0,cv)dio.  (4.7) 

Similarly,  the  zonal-wavenumber  spectrum  of  westward-going  motions  is 

r+co 

rvw(k)  =  /  <  0, io)du>.  (4.8) 

The  spectra  in  the  T /P  data  of  eastward-going  and  westward-going  motions  are  dis¬ 
played  in  figure  4-7.  The  most  striking  feature  in  figure  4-7  is  that  the  energy  level  of 
the  spectrum  of  westward-going  motions  is  higher  than  that  of  eastward-going  ones 
by  a  factor  of  3  at  wavelengths  longer  than  250  km.  At  wavelengths  shorter  than  250 
km,  the  spectra  of  westward-going  and  eastward-going  motions  are  nearly  identical. 

The  general  shape  of  the  zonal- wavenumber  spectrum  of  westward-going  motions 


can  be  written  as 


TrjW(k)  oc  \k\  ^2,  0  <  |&|  <  0.0028  (cycles  per  km) 

oc  \k\~\  0.0028  <  \k\.  (4.9) 

For  the  eastward-going  motions,  it  is 

Tve(k)  oc  |A;|~1//2,  0  <  |A:|  <  0.004  (cycles  per  km) 

oc  |A;|~4,  0.004  <  \k\.  (4.10) 

By  definition,  the  scalar  zonal-wavenumber  spectrum  can  be  written  as 

/+oo  p0  r4*oo 

]{k,uj)duj—  I  uj)duj  +  /  T7?(/c,ct;)(icj  (4.11) 

-oo  J  — oo  J  0 

which  is  symmetrical,  T^(A:)  =  T v(—k),  and  hence  provides  no  information  on  the 
directional  property. 

The  scalar  wavenumber  spectrum  is  the  sum  of  the  wavenumber  spectrum  of 
eastward-going  and  westward-going  motions  in  the  positive  frequencies;  that  is, 

T v(k)  =  T ^(-k)  +  T^A;)  k  >  0.  (4.12) 

The  scalar  zonal- wavenumber  spectrum  is  displayed  in  figure  4-9.  The  spectrum  has 
a  transition  point  around  300  km.  The  spectral  slope  is  about  —1/2  at  wavelengths 
longer  than  300  km  and  increases  to  about  -3  at  wavelengths  shorter  than  300  km, 
which  is  similar  to  that  derived  by  Wunsch  and  Stammer  (1996). 

A  global  wavenumber  spectrum  of  sea  surface  height  has  been  constructed  by 
Wunsch  and  Stammer  (1995).  Subsequently,  Stammer  (1997)  investigated  the  spec¬ 
trum  of  sea  surface  height  region  by  region.  The  global  averaged  wavenumber  spec¬ 
trum  of  sea  surface  height  can  be  written  as 
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Figure  4-6:  The  meridional-wavenumber/frequency  spectra:  (a)  northward-going,  (b)  southward- 
going,  and  (c)  standing  part.  The  annual  cycle  dominates  the  frequency  spectrum  of  stand¬ 
ing  motions. 
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(4,13) 


Tn{k)  oc  |A;|-1/2,  0  <  \k\  <  1/400  (cycles  per  km) 

oc  \k\~5/2,  1/400  <  |A;|  <  1/150  (cycles  per  km) 
cx  |A;|~4,  1/150  <  |A;|. 

The  reader  is  reminded  that  there  is  large  uncertainty  of  the  wavenumber  spectrum  at 
wavelengths  shorter  than  about  100  km.  The  measurements  of  T/P  at  high  wavenum¬ 
bers  are  contaminated  by  measurement  noises  and  aliased  by  small-scale  physical 
processes  such  as  internal  waves  and  internal  tides  (Wunsch  and  Stammer,  1995; 
Stammer,  1997a).  Stammer  (1997a)  found  some  universal  characteristics  for  the  re¬ 
gional  wavenumber  spectral  shape,  which  is  roughly  consistent  with  the  global  one. 
The  energy  level  for  the  sea  surface  height  and  cross-track  velocity  for  wavelengths 
longer  than  about  100  km  increases  greatly  from  the  low  energy  area  in  the  eastern 
and  central  sub-tropical  gyre  toward  the  energetic  western  boundary. 

The  meridional-wavenumber  spectra 

The  meridional-wavenumber  spectra  of  the  northward-going  and  southward-going 
motions  can  be  derived  in  a  similar  way.  Figure  4-8  shows  that  the  meridional- 
wavenumber  spectral  shape  of  northward-going  and  southward-going  motions  are 
similar,  and  the  energy  level  of  southward-going  motions  is  slightly  higher  at  wave¬ 
lengths  longer  than  500  km.  The  scalar  meridional- wavenumber  spectrum  is  displayed 
in  figure  4-10.  It  looks  like  the  scalar  zonal-wavenumber  spectrum  in  figure  4-9. 

The  frequency  spectra 

By  definition,  the  frequency  spectrum  of  eastward-going  energy  can  be  written  as 

Tve(v)  =  f  T  r,{k,u>)dk-,  (4.14) 

Jk>  0 

the  frequency  spectrum  of  westward-going  energy  is 

Tw(a;)  =  [  Trj(k,u)dk]  (4-15) 

Jk<  0 
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and  the  frequency  spectrum  of  standing  waves  is  T1](k  =  0,u>).  The  frequency  spec¬ 
tra  of  the  northward-going,  southward-going  and  standing  motions  can  be  obtained 
similarly. 

The  frequency  spectra  of  eastward-going,  westward-going  and  standing  waves  are 
shown  in  figure  4-11.  The  peak  associated  with  annual  cycle  only  appears  in  the 
spectrum  of  standing  waves.  For  periods  shorter  than  100  days,  the  spectrum  of 
westward-going  waves  is  redder  than  that  of  the  eastward-going  parts.  At  periods 
longer  than  100  days,  the  energy  level  of  westward-going  motions  is  about  three 
times  as  much  as  that  of  eastward-going  motions.  Again,  the  frequency  spectrum  of 
standing  motions  is  dominated  by  seasonal  cycle,  and  the  peak  near  60  days  is  due 
to  the  contamination  by  a  residual  tidal  contribution. 

The  spectral  shape  of  the  frequency  spectra  of  westward-going  motions  can  be 
approximately  written  as 

Tw(w)  oc  l a-1/2,  0  <  u)  <  0.007  (cycles  per  day) 

oc  lu~2,  0.007  <  uj,  (4.16) 

for  the  eastward-going  motions,  it  is 

oc  aA1/2,  0  <  to  <  0.007  (cycles  per  day) 

oc  oT3/2,  0.007  <  u>,  (4.17) 

and  for  the  standing  motions,  it  is 

T^o(w)  oc  ta-3/4,  0  <  to  <  0.005  (cycles  per  day) 

oc  a;-2,  0.005  <  u,  (4.18) 

Figure  4-12  presents  the  frequency  spectra  of  northward- going,  southward-going 
and  standing  waves.  Figure  4-12  is  similar  to  figure  4-11  except  that  the  difference 
between  the  spectra  of  northward-going  and  southward-going  motions  is  much  smaller 
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than  that  between  the  spectra  of  southward-going  and  westward-going  motions. 
The  scalar  frequency  spectrum  is 


f+oo  r+ oo 

/  /  Tv(k,l,uj)dkdl, 

(4.19) 

J — oo  J — OO 

T^cj)  +  Tfj/toG^)  +  T  ^{k  =  0,u>) 

(4.20) 

'Y'ijn  (^0  +  TTt?s(w)  +  Tn(l  =  0,u). 

(4.21) 

The  scalar  frequency  spectrum  is  displayed  in  figure  4-13.  The  global  averaged 
frequency  spectrum  of  sea  surface  height  has  been  constructed  by  Wunsch  and  Stam¬ 
mer  (1995)  and  in  general  the  frequency  spectrum  displayed  in  figure  4-13  is  similar  to 
the  global  averaged  frequency  spectrum  of  sea  surface  height.  On  time  scales  shorter 
than  100  days  the  spectra  approximately  follow  an  lo~2  power  law,  with  an  ta-1/2 
power  law  over  longer  periods  (Wunsch  and  Stammer,  1995). 

Stammer  (1997)  broke  the  global  power  density  spectrum  down  into  regional  struc¬ 
tures.  An  important  result  is  that  a  very  high  degree  of  universality  in  spectral  shape 
was  found — independent  of  the  very  large  variations  in  regional  energy  levels.  Stam¬ 
mer  (1997)  found  three  different  region  types:  (i)  the  energetic  boundary  currents, 
(ii)  the  bulk  of  the  extratropical  basins,  and  (iii)  the  tropical  interior  oceans.  There 
exists  a  pronounced  similarity  in  the  shape  of  all  the  spectra  from  each  dynamical 
category.  In  the  interior  ocean,  the  general  shape  of  the  sea  surface  height  and  sea 
surface  slope  spectra  basically  agrees  with  that  of  the  global  average  but  with  less 
energy.  In  general,  the  slope  spectral  characteristic — with  a  flat  low-frequency  part 
and  a  steeper  decay  at  higher  frequencies-appears  qualitatively  consistent  with  results 
from  moored  current  meter  and  temperature  measurements. 
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west-east  wavenumber  {cycles/km ) 


Figure  4-7:  The  zonal- wavenumber  spectra  of  eastward-going  (solid  line)  and  westward 
going  motions  (dotted  line)  in  the  T/P  data. 


Figure  4-8:  The  meridional-wavenumber  spectra  of  northward-going  (solid  line)  and 
southward-going  motions  (dotted  line). 


Figure  4-11:  The  frequency  spectra  of  eastward-going  (solid  line),  westward-going 
(dotted  line)  and  standing  motions  (dashdot  line).  The  frequency  spectrum  of  stand¬ 
ing  motions  is  dominated  by  annual  cycle. 


Figure  4-12:  The  frequency  spectra  of  northward-going  (solid  line),  southward- going 
(dotted  line)  and  standing  motions  (dashdot  line).  The  frequency  spectrum  of  stand¬ 
ing  motions  is  dominated  by  annual  cycle. 
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seasonal  cycle. 

4.7  Observed  temperature  wavenumber  spectra 

The  data  used  here  is  taken  from  the  repeated  XBT  lines  in  the  North  Pacific  between 
San  Francisco  and  Hawaii.  The  data  have  been  described  in  detail  by  Gilson  et  al. 
(1998).  The  raw  data  were  interpolated  onto  a  grid  with  0.1°  longitude  spacing.  The 
time  mean  of  the  objectively  mapped  data  was  removed  at  each  grid  point.  The  time- 
averaged  wavenumber  spectra  of  the  temperature  perturbations  at  different  depths 
were  then  estimated. 

Figure  4-14  presents  the  time-averaged  wavenumber  spectra  at  different  depths. 
The  general  spectral  shape  is  independent  of  depth  and  can  be  written  as 

T g(k)  oc  |A;|~X/2,  0  <  \k\  <  1/400  (cycles  per  km) 

oc  |*|-5/2,  1/400  <  |*|  <  1/150  (cycles  per  km) 

oc  \k\~\  1/150  <|lfe|,  (4.22) 
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which  is  similar  to  that  of  sea  surface  height  perturbations  (equation  4.13).  In  general, 
the  energy  level  decreases  with  depth.  Roemmich  and  Cornuelle  (1990)  analyzed  XBT 
measurements  in  the  South  Pacific  ocean  and  found  that  the  wavenumber  spectral 
shape  there  is  similar  to  those  displayed  in  figure  4-14. 

Gilson  et  al.  (1998)  studied  the  relationship  of  TOPEX/POSEIDON  altimetric 
height  to  steric  height,  which  is  calculated  from  temperature  and  salinity  measure¬ 
ments  along  a  precisely  repeating  ship  track  in  the  North  Pacific  over  a  period  of  5 
years.  They  found  that  on  wavelengths  longer  than  500  km  altimetric  height  vari¬ 
ability  is  highly  correlated  with  the  variability  of  the  steric  height,  with  a  coherence 
amplitude  of  0.89.  About  65%  of  the  total  variance  is  found  in  the  wavelengths  longer 
than  500  km  in  steric  height. 


Figure  4-14:  The  wavenumber  spectra  of  temperature  averaged  among  depths  (1) 
0-100  m,  (2)  100-200  m,  ...  (8)  700-800  m. 
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4.8  Observed  frequency  spectra  of  velocity  and  tem¬ 


perature 

An  early  discussion  of  frequency  spectra  from  current  meters  was  given  by  Wunsch 
(1981).  Generally  speaking,  such  spectra  of  velocity  show  an  isotropic  high  frequen¬ 
cies  with  a  spectral  slope  of  about  cu-2  followed  by  an  energy-containing  band  towards 
longer  periods  (periods  longer  than  100  days).  In  the  vertical,  the  structure  is  roughly 
consistent  (see  Wunsch,  1997)  with  dominance  of  the  barotropic  and  first  baroclinic 
modes.  Other  results  were  summarized  by  Schmitz  (1978).  These  results  are  modified 
in  proximity  to  major  topographic  features.  Figure  4-14  displays  the  mean  normal¬ 
ized  frequency  spectra  of  horizontal  velocities  in  the  barotropic  and  first  baroclinic 
modes  from  105  current  meter  mooring  measurements,  which  has  been  elaborated 
on  by  Wunsch  (1997).  At  each  site,  the  spectra  have  been  normalized  by  their  total 
energy  before  the  average  is  performed.  The  averaged  frequency  spectra  for  zonal  and 
meridional  velocities  in  the  barotropic  and  first  baroclinic  modes  are  very  similar.  We 
shall  see  later  that  in  most  regions  the  flows  are  isotropic,  in  the  sense  that  the  differ¬ 
ence  between  the  frequency  spectra  of  the  two  components  of  horizontal  velocities  is 
statistically  insignificant.  The  frequency  spectra  of  the  barotropic  and  first  baroclinic 
modes  also  display  similar  structure,  which  implies  that  the  frequency  spectral  shape 
of  horizontal  velocities  is  independent  of  depth.  The  variance-preserving  figures  show 
that  kinetic  energy  are  dominated  by  motions  with  periods  around  100  days. 

Moored  temperature  frequency  spectra  are  not  very  dependent  on  geography  and 
exhibit  a  behavior  similar  to  that  of  velocity.  In  the  regions  away  from  major  topo¬ 
graphic  features,  the  frequency  spectral  shape  of  the  horizontal  velocities  and  tem¬ 
perature  is  independent  of  depth,  and  the  energy  level  of  the  temperature  frequency 
spectra  drops  more  rapidly  with  depth  than  that  of  the  horizontal  velocities. 
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Figure  4-15:  The  mean  normalized  frequency  spectra  of  zonal  (a)  and  meridional 
(b)  velocities  in  the  barotropic  (solid  line)  and  first  baroclinic  modes  (dashed  line) 
from  105  current  meter  measurements  (see  Wunsch  1997).  The  variance-preserving 
form  of  zonal  (c)  and  meridional  (d)  velocities.  At  each  site,  the  spectra  have  been 
normalized  by  their  total  energy  before  the  spatial  average  is  done. 


4.9  Heat  fluxes 


Large  portions  of  the  heat  fluxes  in  the  ocean-atmosphere  system  are  carried  by 
eddies.  The  eddy  heat  fluxes  are  fundamental  to  the  understanding  of  climate  change, 
the  mechanisms  of  eddy  generation,  the  interactions  between  eddy  and  mean  flows 
and  the  parameterization  of  the  eddy  effect  in  ocean  general  circulation  models.  In 
spite  of  the  importance  of  eddy  heat  fluxes,  a  quantitative  global  description  of  eddy 
heat  fluxes  in  the  ocean  is  still  not  available  due  to  the  difficulties  of  observing  the 
ocean  with  sufficient  space  and  time  resolution.  Recently,  Wunsch  (1999a)  estimated 
the  eddy  heat  fluxes  based  on  quasi-global  current  meter  and  temperature  mooring 
records.  He  concluded  that  eddy  heat  fluxes  are  quite  small  in  the  ocean  interior  and 
are  only  significant  near  western  boundary  current  areas.  In  the  western  boundary 
area,  significant  eddy  heat  fluxes  are  confined  to  the  top  1000  m.  Because  there  is  no 
significant  eddy  heat  flux  throughout  much  of  the  ocean  interior,  there  is  no  need  to 
parameterize  the  eddy  heat  flux  there.  The  heat  fluxes  of  the  simple  model  in  Chapter 
3  are  zero.  Therefore  in  the  ocean  interior,  the  data  and  model  are  consistent  to  first 
order. 

Because  of  the  limitation  of  measurements,  the  eddy  heat  fluxes  are  often  esti¬ 
mated  indirectly  through  the  theory  of  baroclinic  instability  or  through  eddy-resolved 
numerical  models.  For  example,  using  four  years  of  sea  surface  height  measurements 
from  the  TOPEX/POSEIDON  spacecraft,  Stammer  (1998)  calculated  the  surface 
eddy  kinetic  energy  and  eddy  scales  under  geostrophic  assumptions,  from  which  the 
global  distribution  of  apparent  eddy  coefficient  was  produced.  Based  on  the  assump¬ 
tion  that  the  eddy  heat  fluxes  are  produced  through  the  baroclinic  instability  of  the 
large  scale  geostrophic  fields,  he  calculated  the  meridional  eddy  heat  fluxes  and  found 
that  strong  meridional  eddy  fluxes  exist  in  the  western  boundary  current  regions  of 
all  oceans  and  in  the  tropics  equatorward  of  about  15°  latitude.  The  significant  eddy 
heat  fluxes  near  the  strong  western  boundary  currents  presumably  result  from  the 
inhomogeneity  of  the  low  frequency  variability  as  described  in  section  4.1. 
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4.10  Anisotropy 


Two  types  of  anisotropy  for  low  frequency  oceanic  variability  are  found  in  this  study. 
First,  the  energy  distribution  is  not  identical  in  all  directions.  For  example,  the  three- 
dimensional  spectrum  of  sea  surface  height  anomaly,  as  described  in  previous  sections, 
shows  that  there  is  more  energy  going  westward  at  large  scales  and  low  frequencies. 
The  second  type  of  anisotropy  is  displayed  in  the  form  of  differences  between  zonal 
and  meridional  components  of  horizontal  velocity,  which  include  the  differences  in 
the  total  energy,  the  differences  in  the  frequency  spectrum  and  the  differences  in 
the  wavenumber  spectrum.  The  difference  between  the  frequency  spectra  of  zonal 
and  meridional  components  of  horizontal  velocity  can  be  found  from  current  meter 
records.  The  wavenumber  spectra  of  zonal  and  meridional  components  of  velocity 
can  be  obtained  only  from  records  of  simultaneous  mooring  deployments. 

The  data  used  here  are  from  the  105  quasi-global  current  meter  moorings,  as 
described  by  Wunsch  (1997).  I  calculated  the  variance  of  zonal  and  meridional  com¬ 
ponents  of  velocity  in  the  barotropic  and  first  baroclinic  mode  for  each  site  and  found 
that  the  variance  of  zonal  and  meridional  velocity  typically  do  not  differ  by  more  than 
a  factor  of  2  in  the  barotropic  and  first  baroclinic  modes.  Because  the  energy  level 
can  not  fully  describe  the  current  meter  records,  I  further  compared  the  frequency 
spectrum  of  zonal  and  meridional  components  in  the  barotropic  and  first  baroclinic 
mode  at  each  site.  It  is  found  that  in  general  the  difference  between  the  frequency 
spectra  of  zonal  and  meridional  components  is  statistically  insignificant  with  respect 
to  a  95%  confidence  interval.  Only  in  the  regions  close  to  boundaries  and  jets,  did 
I  find  significant  anisotropy.  Figures  4-16  and  4-17  are  two  examples  where  signifi¬ 
cant  anisotropy  is  found.  Figure  4-16  shows  that  the  energy  level  of  the  frequency 
spectrum  of  the  zonal  component  is  over  one  order  of  magnitude  higher  than  that 
of  the  meridional  component  at  periods  longer  than  100  days  in  both  the  barotropic 
and  first  baroclinic  modes.  The  result  in  figure  4-17  is  the  opposite;  at  periods  longer 
than  50  days,  the  energy  level  of  the  meridional  component  is  significantly  higher 
than  that  of  the  zonal  component  in  the  two  modes.  The  data  used  in  figures  4-16 
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and  4-17  are  taken  inside  the  Gulf  Stream  and  near  the  eastern  boundary  of  the  North 
Pacific,  respectively.  Therefore,  the  difference  between  the  frequency  spectrum  of  the 
zonal  and  meridional  components  in  figure  4-16  is  presumably  related  to  the  strong 
western  boundary  current,  the  Gulf  Stream.  The  difference  in  figure  4-17  might  be 
due  to  the  meridional  California  Current  system.  The  results  in  figure  4-16  and  4-17 
suggest  that  the  dynamics  in  the  region  close  to  jets  and  boundaries  are  distinct.  The 
reason  for  the  difference  between  the  zonal  and  meridional  components  as  displayed 
here  awaits  further  study. 
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Figure  4-16:  Frequency  spectra  for  the  zonal  (solid  line)  and  meridional  (dotted 
line)  component  of  horizontal  velocity  in  the  barotropic  and  first  baroclinic  mode  at 
38.8°N,  68.1°W.  The  spectra  are  estimated  using  D.  Thompson’s  multitaper  method 
(Percival  and  Walden,  1993).  The  ratio  of  the  variance  of  zonal  component  to  that 
of  meridional  component  in  the  barotropic  and  first  baroclinic  mode  is  3.5  and  2.7, 
respectively. 
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Figure  4-17:  Same  as  figure  4-16  except  at  38.2°N,  124.4°W.  The  ratio  of  the  vari¬ 
ance  of  zonal  component  to  that  of  meridional  component  in  the  barotropic  and  first 
baroclinic  mode  is  0.34  and  0.25,  respectively. 
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Chapter  5 


Energy  distribution  in  k,  l  and  u 
space  part  I:  a  scalar  form 

5.1  Fitting  <!>(&,  /,  tu,  n ,  </>,  A)  from  observations 

Based  upon  the  summary  in  chapter  4  of  observed  frequency  and  wavenumber  spectra, 
I  will  choose  forms  for  the  separated  functions  making  up  the  full  $(k,l,u,n,(j>,  A). 
As  demonstrated  in  chapter  3,  once  the  three-dimensional  spectrum  for  pressure  of 
each  mode  is  known,  one  can  estimate  A)  through  the  equation  (3.8). 

Unfortunately,  none  of  the  full  three-dimensional  spectrum  of  each  mode  for  any  vari¬ 
able  is  available.  Because  the  observed  one-dimensional  frequency  and  wavenumber 
spectral  shapes  are  approximately  independent  of  depth,  one  can  infer  that  the  fre¬ 
quency  and  wavenumber  spectral  shapes  of  each  mode  are  approximately  identical, 
and  thus  the  formula  for  $>(k,  l,LU,n,  (f>,  X)  will  be  the  same  for  each  vertical  mode. 
Because  the  full  three-dimensional  spectrum  is  too  complicated  to  be  presented  in 
an  analytical  form,  I  start  here  with  the  simple  form  as  given  by  equation  (4.3).  It 
is  assumed  that  the  frequency  spectral  shape,  wavenumber  spectral  shape,  vertical 
structure  and  energy  level  are  separable.  The  separated  form  in  equation  (4.3)  is  in¬ 
capable  of  differentiating  the  westward-going  energy  from  the  eastward-going  energy. 
A  more  sophisticated  frequency  and  wavenumber  coupled  form  will  be  produced  in 
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chapter  7.  In  the  following,  I  will  choose  the  form  of  Bn(k),  Cn(l ),  Dn(u)),E0(n )  and 
I ((f),  A)  in  equation  (4.3)  according  to  the  observations.  Because  most  of  the  region 
is  dominated  by  the  first  few  modes  as  seen  in  the  last  chapter,  only  barotropic  and 
first  two  baroclinic  modes  are  included  in  my  model:  n=0,  1,  2.  The  form  of  Bn(k ), 
Cn(l)  and  Dn(u) )  are  chosen  to  fit  the  globally  averaged  wavenumber  and  frequency 
spectral  shape,  respectively. 

Choosing  the  Bn(k )  to  fit  the  observed  zonal-wavenumber  spectrum 

The  Bn(k )  is  chosen  according  to  the  globally  averaged  wavenumber  spectral  shape 
of  sea  surface  height  from  T /P  measurements: 


a0\k\2 


Bn(k)  =  { 


ai\k\  ly/2 

a2\k\~5/2 


if  0  <  \k\  <  k0 
if  ko  <  |&|  <  ki 
if  ki  <  |&|  <  k2 
if  \k\  >  k2  , 


(5.1) 


where  k2  =  0.008  and  k\  =  0.0025  (cycles  per  km).  For  continuity  of  Bn(k)  at  k2 
and  k\,  a2  —  k^^2  and  aq  =  k^^k^2.  In  the  ocean,  the  longest  wavelength  is 
limited  by  the  width  of  the  ocean  basin  or  the  perimeter  of  the  earth  where  the 
ocean  circulates  around  the  earth  without  barriers.  Therefore,  there  is  no  energy 
at  wavenumbers  0  <  \k\  <  ko,  where  ko  corresponds  to  the  wavenumber  of  longest 
wave  in  the  ocean.  If  the  spectrum  at  wavenumbers  0  <  |&|  <  k0  is  chosen  to  be 
zero,  it  will  be  discontinuous  at  |fc|  =  ko,  and  thus  the  so-called  Gibbs  phenomenon 
will  arise.  In  order  to  moderate  the  Gibbs  phenomenon,  the  spectrum  is  chosen  to 
diminish  smoothly  to  zero  from  k  =  k0  to  k  =  0,  while  it  ensures  that  the  energy  at 
wavelengths  longer  than  the  width  of  the  ocean  basin  is  negligible.  As  displayed  in 
equation  (5.1),  the  wavenumber  spectrum  is  simply  chosen  to  be  proportional  to  \k\2 
at  wavenumbers  0  <  \k\  <  k0.  For  continuity  of  Bn(k)  at  k0,  «0  =  kfi^^k^kfifi^2 .  For 
simplicity,  ko  is  assumed  to  be  0.0001  cycles  per  km.  The  corresponding  wavelength 
of  k0  is  10,000  km,  which  is  about  as  long  as  the  width  of  the  North  Pacific  ocean. 
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Equation  (5.1)  gives 


/+OO 

Bn(k)dk 

-OO 


5.1  x  107. 


(5.2) 


Choosing  the  Cn(l)  to  fit  the  observed  meridional-wavenumber  spectrum 
The  wavenumber  spectra  of  sea  surface  height  from  T/P  measurements  suggest 
that  Cn(l )  has  the  same  form  as  Bn{k). 

Choosing  the  Dn(oj)  to  fit  the  observed  frequency  spectrum 

The  Dn(w)  is  chosen  according  to  the  globally  averaged  frequency  spectral  shape 
of  horizontal  velocity: 


f  7oM~1/2 


U) 


if  0  <C  \uj  |  ^  cvq 
if  \lo\  >  u0  , 


(5.3) 


where  u;0  =  0.01  (cycles  per  day).  For  continuity  of  Dn{k )  at  cv0,  q0  =  co0  3//2.  Equation 
(5.3)  gives 


/+oo 

Dn(u))du>  =  600. 

-OO 


(5.4) 


Now  that  the  wavenumber  and  frequency  spectral  shapes  have  been  chosen,  the 
variance  for  each  mode  can  be  obtained 


/+oo  r-t oo  r+oo 

/  /  $(k,fiLO,n,(f),  \)dkdldoj 

-oo  J  —  oo  J  —  oo 

=  X)E0(n)  y+^Bn(k)d^  [iyn(l)dl}  [/_+< 

=  1.6  x  1O187(0,  A)Eo(rr). 

The  energy  level  of  each  mode 

The  total  kinetic  energy  per  unit  surface  area  for  each  mode  is 


(5.5) 


/+ oo  r+oc  /*- f  oo 

/  /  27T2(k2  +  l2)$>(k,l,uj,n,<i),\)dkdlduj 

-oo  J ~ oo  J — oo 
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f  /*+ oo  If  f+oo  1  r  r+oo 

=  2n2I((t>,\)E0{n)  /  k2Bn(k)dk  /  Cn(l)dl  /  £>„(u;)d< 

L*/  — oo  J  U-oo  J  .J —oo 

f  Z’+OO  ]  f  r+oo  If  /»+oo 

+2tv2 I ((f),  X)E0(n)  /  £«(*)<**  /  l2Cn{l)dl  /  Dn(w)da; 

!_«/ —  oo  J  L«/— oo  J  [_</ —  oo 

=  6.3  x  lO14/(0,A)£o(n).  (5.6) 


As  seen  in  section  4.3,  the  vertical  structure  of  eddy  kinetic  energy  can  be  simply 
represented  by  the  barotropic  and  first  two  baroclinic  modes  in  the  middle  of  the 
ocean,  and  roughly  speaking,  the  kinetic  energy  of  the  first  three  modes  is  in  the  ratio 
of  1 :1 : 1/2.  Accordingly,  if  we  choose  £o(0)  =  1,  then  £0(1)  =  1  and  E0{2)  =  0.5. 
Choose  I{<p,  A)  to  fit  the  observed  surface  kinetic  energy 
The  surface  kinetic  energy  for  the  model  is 


f  r+ OO  /*+ oo  r+oo 

Ek((f>,X,z  =  0)  =  J2F2(z  =  0)  /  /  /  2ir2(k2  +  l2)$(k,l,u,n,(l),\)dkdldi 

„_n  U-oo  J  —  oo  J— oo 

(5.7) 


n=0 

n=2 


E  FKZ  =  0 )KE(n)  =  2.1  x  10 12I(<f>,  A). 

n=0 


Equation  (5.7)  yields 


I(<f>,  A) 


Ek(<f>,  \z  —  0) 
2.1  x  1012  ' 


(5.8) 


The  empirical  formulas  for  the  surface  kinetic  energy  in  the  North  Pacific  and  the 
North  Atlantic  are  given  by  equations  (4.1)  and  (4.2),  respectively. 

The  potential  energy  per  unit  surface  area  for  each  mode 


1  f+ oo  /*+ oo  r-foo 

PE(n,<j>,  A)  =  — —  /  /  /  $(k,l,w,n,(f>,\)dkdldcu 

ZItt  J— 00  J —oo  J —oo 


7.8  x  lO17/(0,A)£o(n) 


R2 

±Ln 


(5.9) 


The  ratio  of  kinetic  energy  to  potential  energy  per  unit  surface  area  for  first  baro¬ 
clinic  mode 


KE{  1,&  A) 
A) 


=  8.1  x  10“4^. 


(5.10) 
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If  we  take  /  =  10  4s  2,  which  is  a  typical  value  of  /  at  mid-latitude,  equation  (2.21) 
gives  Ri  =  25km.  Then 


KE{  UA) 
PE(  1,6  A) 


=  0.51. 


(5.11) 


77m  zonal-wavenumber  spectrum  of  sea  surface  height 


T„(M,A)  =  £ 


^  f2F%(z  =  0)  /-+00  /-+00 
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0  f  r+oo  ]  r  /*+ oo  1 

E  *?(*  =  0)^o(n)  /  Cn{l)dl  /  7}„(u;)du;  7?n(7).  (5.12) 
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r-t oo  r-Hoo 

Mote  the  consistency  of  units  for  different  variables.  /  /  $(k,l,u>,n,  (j),  \)dldu) 

has  the  units  of  (cm2/s2)m/(cycles/km)3,  the  units  of  F2(z)  are  m_1,  the  units  of  / 
are  s-1,  and  the  units  of  g  are  km  s-2.  If  we  take  /  =  2Q.sin<t)  =  1.45  x  lO_4sm0  s_1 
and  j  =  9.8x  10_3km  s-2,  then 


^r)(k,(j>,  A) 


/2/(0,  A)  x  5.1  x  107  x  600 


[F02(z  =  0)  +  F?{z  =  0)  +  0.5  x  Fpz  =  O^B^k) 


—  2.2  x  104 sin2 <f  7(0,  X)B\ (k)  [cm2/(cycles/km)]. 


(5.13) 


The  meridional-wavenumber  spectrum  of  sea  surface  height 
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=  t-Mn.  2  Fl(z  =  0)B„(n)  [£*(*)*]  £°°n„(w)da>]  C„(() 
—  2.2  x  104sm2^)  7(0,  A)Ci(/)  [cm2/(cycles/km)j.  (5-14) 


The  variance  of  sea  surface  height 
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=  1.1  x  1012smV  7(0,  A) 


=  0.52 sin2f>  Ek{(j),  A,  2  =  0)  [cm2  . 


(5.15) 


The  zonal-wavenumber  spectrum  of  temperature 
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=  6.4  x  10Asin2f)I(f),X)[G\(z)  +  0.5  x  G\(z)](-^)2  Bi(k)  [°C2/(cycles/km)],(5.16) 

where  the  units  of  G2(z)  are  s4m~3  and  the  units  of  dOo/dz  are  °C/m. 

The  meridional-wavenumber  spectrum  of  temperature 
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rjfj^  r  r+00  1  r  r+oo 

=  f2I(<t>,\)(?P-)2[G2(z)  +0.5  x  G2(z)]\  Bx(k)dk  /  Dx(u)du  Cx(l ) 

OZ  U —00  J  U— 00 

30 

=  6.4  x  104sm2^  1(f),  A) [G\ (z)  4-  0.5  x  G\(z)}(-^-)2Cx(l)  [°C2/(cycles/km)],(5.17) 

. .  oz 

where  the  units  of  G2(z)  are  s4m~3  and  the  units  of  dO^/dz  are  °C/m. 

The  frequency  spectrum  of  sea  surface  height 


T,(w,0,  A) 


n— 2  —  01  r+00  r+00 

V  - — Z-  I  /  $(k,l,u,n,<f>,\)dhdl 

n= 0  9  J~°°  J~°° 

— •  E  S2(z  =  O)So(n)  [£~BB(*)d*]  |^CB(Q<fl]  A,(w) 


=  1.9  x  109sm2 <?!>/(</>,  A)-Di(u;)  [cm2/ (cycles/km)]  .  (5.18) 

The  frequency  spectrum  of  zonal  velocity 

Tt~^  f-00  r+oo 

Tu(w,z,(j),\)  =  /  £n2l2F2(z)$(k,l,u,n,f),  X)dkdl 

71  =  0  •/-°° 

=  Ait2 1(f>,  A)  E0(n)F2(z  =  0)  [|_+Jsn(fc)dA;]  [|_J2Cn(0^]  A»(w) 

=  1.1  x  1012  [F02(z)  +  F2(z)  +  0.5  x  S|(z)J  Di(u>)  [cm2/s2/(cycles/day)j,  (5.19) 
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where  the  units  of  F2(z)  are  m-1. 

The  frequency  spectrum  of  meridional  velocity 

n=2  n+oc  r- foo 

T v{aj,z,(j>,\)  =  J2  /  4ir2k2F2(z)$(k,  l,u,n,  0,  X)dkdl 

n=0J~°° 

n—2  r  />-foo  r  /»+ oo 

=  4iT2I((j),  A)  F2(z)E0(n)  /  k2Bn(k)dk  /  Cn(l)dl  AM 

n=0  U-oo  J  U-oo  J 

=  1.1  x  1012  jF02(z)  +  F2(z)  +  0.5  x  F2(z)]  AM  [cm2/s2/(cycles/day)],  (5.20) 

where  the  units  of  F2(z)  are  m_1.  Note  that  the  frequency  spectra  of  zonal  and 
meridional  velocities  are  the  same. 

The  frequency  spectrum  of  temperature 

F&(oj,  z,  (f>.  A)  =  /2G'2(z)(— )2  /  /  <b(k,l,u,n,(f>,  \)dkdl 

n~Q  uZ  J  —oo  J  —  oo 

=  /2/(</>,  A)(^)2[G2(z)  +  0.5  x  G2(z)j  (*)<**]  [|+°°C1(0d/|  AM 

=  5.4  x  lO9sm20  /AA)[G?M  +  0.5  x  G2(z)](^)2AM  [°C2/(cycles/day)]{5.21) 

where  the  units  of  G2(z)  are  s4m-3,  the  units  of  d90/dz  are  °C/m. 

5.2  Model  and  data  comparison 

In  the  above  section,  I  have  specified  the  formula  for  each  component  of  <!>(&,  I,  u>,  n,  <f>,  A) 
Bn(k),  Cn(l),  Dn(u>),  E0(n )  and  I(<p,  A),  based  on  different  observations.  There  may 
be,  however,  some  incompatibilities  among  different  observations.  For  example,  the 
vertical  structure  of  the  model  is  specified  according  to  the  observed  vertical  structure 
of  kinetic  energy.  However,  it  is  unclear  whether  the  vertical  structure  of  the  temper¬ 
ature  (potential  energy)  is  dynamically  consistent  with  that  of  the  kinetic  energy.  An¬ 
other  example  is  the  wavenumber  spectral  shape.  Both  T/P  and  XBT  measurements 
can  provide  estimates  of  wavenumber  spectral  shape.  If  they  are  not  dynamically 
consistent,  the  function  $(k,  l,  n,  co,  A)  through  which  diverse  observations  can  be  rec- 


onciled  will  not  exist.  In  this  section,  I  will  make  model/data  comparisons  to  see  how 
well  each  model  spectrum  fits  the  corresponding  observation  with  the  $(&,  l,  oj,  n,  0,  A) 
we  chose.  We  do  not  expect  the  model  to  fit  the  observations  perfectly,  because  the 
different  wavenumber,  frequency  and  vertical  mode  structure  found  in  different  places 
have  been  neglected,  the  seasonal  cycle  is  not  considered  either,  etc. 

Zonal- wavenumber,  meridional- wavenumber  and  frequency  spectrum  of 
sea  surface  height 

The  model  and  observed  zonal- wavenumber,  meridional- wavenumber  and  frequency 
spectra  are  plotted  in  figure  5-1,  5-2  and  5-3.  The  observed  spectra  come  from  the 
area  of  25°  -  50°N  and  195°  -  225°E  with  a  center  latitude  at  37.5°N,  which  have  been 
described  in  detail  in  section  4.4.  The  energy  level  inside  this  area  is  relatively  uni¬ 
form.  According  to  the  equation  (4.1),  the  average  value  Ek(<f>,  A,  z  =  0)  in  this  area 
is  about  94  cm2/s2;  therefore,  I((f), A)  =  4.5  x  10-n.  Figures  5-1,  5-2  and  5-3  show 
that  the  model  spectra  are  reasonable  fits  to  the  observations  in  general.  One  major 
difference  between  the  observations  and  the  model  is  that  at  wavelengths  shorter  than 
200  km,  the  observed  wavenumber  spectrum  is  redder  than  the  model  spectrum.  As 
discussed  in  section  4.4,  there  is  some  uncertainty  of  the  observed  wavenumber  spec¬ 
tral  shape  at  wavelengths  less  than  200  km.  Another  obvious  difference  is  that  the 
observed  frequency  spectrum  is  slightly  redder  than  the  model  spectrum  at  periods 
longer  than  100  days. 

Temperature  wavenumber  spectrum 

The  observed  temperature  wavenumber  spectrum  in  figure  5-4  is  from  the  XBT 
data  in  the  eastern  North  Pacific  described  in  section  4.7.  The  XBT  line  ranges  from 
(236°E,  38°N)  to  (202°E,  21°N)  and  the  center  latitude  is  30°N.  The  surface  eddy 
kinetic  energy  along  this  XBT  line  is  relatively  homogeneous.  According  to  equa¬ 
tion  (4.1),  the  mean  value  of  surface  kinetic  energy  along  this  XBT  line  is  about  150 
cm2/s2.  Because  the  vertical  gradient  of  mean  temperature  changes  strongly  with 
location  in  the  upper  ocean,  which  violates  the  assumption  of  horizontal  homogene¬ 
ity,  we  compare  the  model  spectrum  with  the  observations  at  the  depth  of  700  m  as 
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a  representative  result.  The  vertical  gradient  of  mean  temperature  along  the  XBT 
line  at  this  depth  changes  relatively  less  sharply.  The  mean  value  of  d90/dz  is  about 
0.0041°C/m.  With  this  choice  of  energy  level  and  vertical  gradient  of  mean  tem¬ 
perature,  the  model  temperature  wavenumber  spectrum  is  plotted  in  figure  5-4.  As 
displayed  in  figure  5-4,  the  model  spectrum  fits  the  observations  very  well  in  terms  of 
spectral  shape  and  energy  level. 

Frequency  spectrum  of  horizontal  velocities  and  temperature 

One  set  of  current  meter  mooring  data  is  analyzed  here,  which  is  located  at 
(39.5°N,  232. 3°E)  inside  the  low  energy  area  in  the  eastern  North  Pacific.  It  has 
current  and  temperature  meters  at  nominal  depths  of  160,  580,  595,  1230,  3000  and 
3800  m.  The  measurements  at  3800  m  are  neglected  because  the  signal  for  tem¬ 
perature  perturbations  is  very  weak  and  the  record  is  primarily  dominated  by  the 
noise.  From  equation  (4.1),  Ek(39.5°N,232.3°E,  z  =  0)  =  86cm2/s2.  Therefore, 
7(39. 5°iV,  232° E)  =  4.1  x  10~u.  The  vertical  gradients  of  the  mean  temperature 
at  depths  of  160,  580,  595,  1230, and  3000  are  0.020,  0.0042,  0.0040,  0.0019  and 
0.00024°C/m,  respectively.  The  model  and  observed  spectra  are  plotted  in  figure  5-5. 
The  general  agreement  between  the  model  and  the  observed  spectra  is  quite  pleasing. 

Ratio  of  kinetic  energy  to  potential  energy 

The  ratio  of  kinetic  energy  to  potential  energy  of  the  first  baroclinic  mode  is  0.51 
in  mid-latitudes  (equation  5.11).  This  value  is  consistent  with  the  result  obtained  by 
Wunsch(1999b)  from  current  meter  and  temperature  mooring  measurements. 

5.3  Summary 

In  this  chapter,  a  regional  model  for  the  energy  density  $(k,l,u j,n,<f>,X)  of  the  low 
frequency  variability  in  the  ocean  is  produced.  Only  the  first  three  modes  are  included 
in  my  model:  the  barotropic  and  first  two  baroclinic  modes.  For  each  mode  I  choose 

A)  =  Bn{k)Cn{l)Dn{(jj)E0{n)I{<j),X),  (5.22) 
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and  based  on  various  observations,  I  obtain: 

(1)  the  meridional-  and  zonal-wavenumber  spectral  shape 


Cn(l)  =  Bn(k)  = 


2.2  x  1018|A;|2  if  0  <  |£;|  <  k0 
2.2  x  108| A:|— if  k0  <  \k\  <  k\ 
1.4  x  103|&;|-5/2  if  k\  <  \k\  <  k2 
|/c|-4  if  \k\  >  k2  , 


where  k2  =  0.008,  kx  =  0.0025  and  k0  =  0.0001  (cycles  per  km). 

(2)  the  frequency  spectral  shape 


Dn{oj)  = 


1000|u;|  */2  if  0  <  |oj|  <  u 0 


if  u>\  >  ujq  , 


where  u0  —  0.01  (cycles  per  day). 

(3)  the  partition  among  vertical  modes 


(5.23) 


(5.24) 


E0{0)  =  1,E0(1)  =  1,E0(2)=0.5. 


(5.25) 


(4)  the  energy  magnitude 


m  a) 


Ekifa  A,z  =  0) 
2.1  x  1012 


(5.26) 


The  empirical  formula  for  the  surface  kinetic  energy  in  the  North  Pacific  is 

=°)  -  30+il^+1000exp{-[iAw£+Mir£l} 

(5'27) 

where  A  is  the  longitude  ,  <f>  is  the  latitude  and  the  units  of  Ek  are  cm2/s2. 
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The  empirical  formula  for  the  surface  kinetic  energy  in  the  North  Atlantic  is 


Ek{(j>,  A,  z  =  0)  =  50  T 


!L  +  1„ooexpHAri25)!+(izi!)!|}_ 


280expH<A^  +  (* 


sin  (f> 

16)2 


400 


80 


2000 


200 


B-ieM-[M  +  M1}  (5,8) 


900 


50 


where  the  units  of  E*  are  cm2/s2. 

The  frequency  spectral  shape,  wavenumber  spectral  shape  and  vertical  structure 
are  universal  in  my  model;  only  the  energy  level  is  a  function  of  geography.  The 
parameters  for  the  wavenumber  and  frequency  spectral  shape  in  my  model  are  pri¬ 
marily  chosen  to  fit  the  globally  averaged  spectra  of  the  observations.  This  empirical 
model  is  a  good  first-order  approximation  to  the  true  ocean,  because  the  observed 
frequency,  wavenumber  and  mode  structure  displays  strikingly  universal  characteris¬ 
tics.  Nonetheless,  we  need  to  keep  in  mind  that  the  frequency,  wavenumber  and  mode 
structure  do  show  some  variations  from  region  to  region,  even  if  these  variations  are 
very  small  compared  with  the  variations  of  eddy  kinetic  energy  with  geography.  These 
changes  of  frequency,  wavenumber  and  vertical  structure  with  location  are  of  great 
importance  to  the  understanding  of  eddy  generation  mechanism  (Stammer,  1997a). 
Based  on  different  observations,  a  rough  estimate  of  the  uncertainty  of  my  empirical 
model  is  provided  here.  First,  the  transitional  wavenumber  k\  in  equation  (5.1)  is 
assumed  to  be  universal,  and  ki  =  0.0025  cycles/km.  The  study  by  Stammer  (1997a) 
shows  that  k\  is  a  function  of  latitude,  and  it  increases  with  latitude,  ranging  from 
0.002  cycles/km  at  15°N  to  0.005  cycles/km  at  60°N.  Let  us  assume  Bn{k)  oc  kp  and 
Dn(u)  oc  (jf .  The  values  of  p  and  r  at  different  wavenumber  and  frequency  bands  in 
my  model  are  given  in  equation  (5.1)  and  (5.3),  respectively.  In  general,  the  uncer¬ 
tainty  of  p  is  greater  at  high  wavenumbers.  At  wavelengths  longer  than  about  400  km, 
the  global  average  value  of  p  is  -0.5.  The  study  by  Stammer  (1997a)  showed  that  the 
spectral  slope  p  at  this  wavenumber  band  is  nearly  independent  of  location  (his  figure 
16),  and  that  p  changes  from  about  —3.0  to  —2.2  over  the  wavelengths  between  about 
100  and  400  km  from  region  to  region.  There  is  great  uncertainty  associated  with  the 
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wavenumber  spectral  shape  from  T /P  measurements  at  wavelengths  less  than  about 
100  km  because  of  the  measurement  noise  and  aliasing  of  interval  waves.  Wunsch 
and  Stammer  (1995)  obtained  a  fc~4//3  power  law  at  this  high  wavenumber  band  and 
argued  that  this  is  dominated  by  noise.  The  wavenumber  spectra  of  high-quality 
XBT  measurements  in  chapter  4  show  that  the  slope  is  about  —4  at  wavelengths 
less  than  150  km.  Overall,  the  change  of  p  with  location  is  within  about  15%  in 
respect  to  the  global  averaged  value  of  p  except  for  the  wavelengths  less  than  100 
km.  There  is  some  uncertainty  of  the  frequency  spectral  shape  also.  The  globally 
averaged  spectral  slope  of  horizontal  velocity  is  —0.5  at  periods  longer  than  100  days; 
however,  the  spectral  slope  of  sea  surface  hight  is  about  —0.8  and  some  local  current 
meter  measurements  show  that  the  spectrum  of  horizontal  velocity  is  nearly  white  at 
periods  longer  than  100  days.  Although  the  study  by  Wunsch  (1997)  showed  that  on 
average  the  partitioning  of  kinetic  energy  among  barotropic  and  first  baroclinic  mode 
is  roughly  equal,  the  exact  value  for  it  is  a  function  of  location.  For  example,  the 
mid- Atlantic  ridge  region  shows  strong  first  baroclinic  mode  contribution.  Interested 
readers  should  consult  Wunsch  (1997). 
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west-east  wavenumber  (cycles/km) 

Figure  5-1:  Zonal-wavenumber  spectrum  of  sea  surface  height.  Solid  line  is  the  ob¬ 
served  spectrum.  Dotted  line  is  the  corresponding  model  spectrum. 


north-south  wavenumber  (cycles/km) 


Figure  5-2:  Meridional-wavenumber  spectrum  of  sea  surface  height.  Solid  line  is  the 
observed  spectrum.  Dotted  line  is  the  corresponding  model  spectrum. 


frequency  (cycles/day) 


Figure  5-3:  Frequency  spectrum  of  sea  surface  height.  Solid  line  is  the  observed 
spectrum.  Dotted  line  is  the  corresponding  model  spectrum.  There  is  a  striking  peak 
at  the  annual  period  on  the  observed  spectrum. 


Figure  5-4:  Zonal-wavenumber  spectrum  of  temperature.  Solid  line  is  the  observed 
spectrum.  Dotted  line  is  the  corresponding  model  spectrum. 


(cm/s)  /CYCLE/DAY 
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CYCLES/DAY 
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(v) 

Figure  5-5:  Frequency  spectra  of  zonal  (a),  meridional  velocity  (b)  and  temperature 
(c)  at  depths  of  160  m  (i),  580  m  (ii),  595  m  (iii),  1230  m  (iv)  and  3000  m  (v).  Solid 
line  is  the  observed  spectrum,  and  the  dotted  line  is  the  model  spectrum. 


Chapter  6 


Application 


A  first  order  spectral  description  of  low  frequency  oceanic  variability  has  been  con¬ 
structed;  however,  I  have  arbitrarily  ignored  the  contributions  of  higher  modes, 
the  different  structure  found  in  intense  current  systems,  the  equatorial  region,  etc. 
Nonetheless,  the  present  first  approximation  does  have  some  useful  applications.  The 
followings  describe  four  examples  of  the  uses  of  the  spectra. 


6.1  Covariance  functions 

The  Fourier  transform  of  the  wavenumber  (frequency)  spectrum  is  the  spatial  (tem¬ 
poral)  covariance.  For  example,  the  zonal  autocovariance  of  temperature  is 

n=2  .-|-oo  r+oo  /-+oo 

Re{rx,(j>,\,z)  =  Y,  /  /  I  *«(M,  u ;,  z,  n)|2<h(A;,  l ,  u>,  n,  <j>,  \)cos(2irkrx)dkdldu. 

n_0  J  —  OO  J  —  OO  J  —  oo 

(6.1) 

Knowledge  of  the  covariance  function  is  useful  for  many  purposes,  such  as  objective 
mapping  and  the  estimation  of  integral  scales  of  variability.  In  order  to  conveniently 
compare  data  with  different  scales,  the  autocorrelation  function  is  often  calculated  in¬ 
stead.  The  autocorrelation  of  temperature  as  a  function  of  zonal  separation  distance 
is  displayed  in  figure  6-1.  The  temperature  autocorrelation  function  shows  an  expo¬ 
nential  decay  with  e-folding  spatial  scales  of  about  120  km.  Because  the  meridional 
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wavenumber  spectrum  of  temperature  is  the  same  as  the  zonal  wavenumber  spectrum, 
the  meridional  autocorrelation  function  and  the  zonal  autocorrelation  function  are  the 
same  also.  Figure  6-2  shows  the  autocorrelation  of  zonal  and  meridional  velocity  as 
a  function  of  zonal  separation  distance.  The  autocorrelation  function  of  meridional 
velocity  has  a  zero  crossing  at  50  km.  The  most  striking  feature  in  figure  6-2  is  that 
the  autocorrelation  function  of  zonal  velocity  decays  with  zonal  lag  much  faster  than 
that  of  meridional  velocity,  which  means  that  along  zonal  direction,  zonal  velocity 
has  longer  spatial  scales  than  meridional  velocity.  This  is  due  to  the  fact  that  zonal 
wavenumber  spectrum  of  zonal  velocity  is  redder  than  that  of  meridional  velocity,  the 
ratio  of  zonal  wavenumber  spectrum  of  zonal  velocity  to  that  of  meridional  velocity 
is  proportional  to  k~2.  The  meridional  autocorrelation  functions  are  the  same  as  the 
zonal  autocorrelation  functions  except  that  autocorrelation  function  of  zonal  veloc¬ 
ity  and  autocorrelation  function  of  meridional  velocity  are  interchanged.  My  result  is 
consistent  with  the  theory  of  turbulence.  Batchelor  (1953)  found  that  for  incompress¬ 
ible  fluid,  the  longitudinal  velocity  autocorrelation,  which  stands  for  autocorrelation 
of  the  velocity  parallel  to  the  separation  vector,  is  always  higher  than  the  lateral 
velocity  autocorrelation  which  represents  autocorrelation  of  the  velocity  normal  to 
the  separation  vector.  As  discussed  at  the  end  of  chapter  5,  the  universal  wavenum¬ 
ber  spectral  shape  in  my  model  has  some  uncertainty.  Overall,  the  uncertainty  of 
wavenumber  spectral  slope  is  within  15%  in  respect  to  the  universal  spectral  slope  in 
my  model  at  wavelengths  longer  than  100  km.  Figure  6-3  shows  how  sensitive  the 
autocorrelation  function  is  to  the  wavenumber  spectral  shape.  As  displayed  there, 
the  variations  of  the  autocorrelation  associated  with  the  changes  of  spectral  shape 
are  less  than  0.1  everywhere.  The  most  sensitive  place  is  at  the  zonal  lag  between 
200  and  400  km  where  the  maximum  variations  occur.  The  reader  is  reminded  that 
the  energy  of  the  temperature  is  concentrated  at  the  spatial  scale  around  200  to  400 
km.  The  dashdot  line  in  figure  6-3  has  higher  autocorrelation  for  the  same  separation 
distance,  because  it  corresponds  to  the  redder  wavenumber  spectrum. 
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Figure  6-1:  Temperature  autocorrelation  as  a  function  of  zonal  separation  distance. 


Figure  6-2:  Autocorrelation  of  zonal  velocity  (solid  line)  and  meridional  velocity 
(dashdot  line)  as  a  function  of  zonal  separation  distance. 


Figure  6-3:  Uncertainty  of  the  autocorrelation  function.  The  solid  line  represents 
the  autocorrelation  function  corresponding  to  the  universal  wavenumber  spectrum  in 
my  model.  The  dashdot  and  dotted  line  correspond  to  the  wavenumber  spectrum 
with  a  slope  of  15%  more  and  15%  less  than  the  slope  of  my  universal  spectrum, 
respectively. 

6.2  Variability  of  volume  flux  and  heat  content 

Large-scale  long-term  change  of  temperature  in  the  ocean  has  attracted  much  atten¬ 
tion  in  the  last  decade.  The  hydrographic  section  along  24°N  in  the  North  Atlantic 
provides  unique  data  to  detect  climate  change  of  large  scale  temperature  and  volume 
transport  in  the  ocean,  because  it  has  been  repeated  three  times  during  1957,  1981 
and  1992.  By  comparing  the  1981  and  1957  surveys,  Roemmich  and  Wunsch  (1984) 
first  identified  significant  warming  in  an  ocean-wide  band  from  700  to  3000  m  with  a 
maximum  temperature  difference  of  0.2°C.  Recently  Parrilla  et  al.  (1994)  and  Arbic 
and  Owens  (2000)  used  all  three  surveys  along  24°N  and  found  that  the  waters  be¬ 
tween  800  and  2500  m  have  consistently  warmed  over  the  past  35  years.  In  climate 
study,  the  natural  variability  is  often  regarded  as  noise  and  the  long  term  trend  is 
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considered  to  be  the  signal.  It  is  highly  desirable  not  only  to  estimate  the  difference 
between  measurements  taken  at  two  different  times,  but  also  to  determine  whether 
such  a  difference  is  significant  relative  to  natural  variability. 

The  volume  transport  of  the  ocean  is  often  estimated  from  hydrographic  mea¬ 
surements  through  inverse  modeling  (Wunsch,  1996).  It  varies  with  time,  and  the 
variability  depends  on  the  horizontal  distance  of  the  hydrographic  section.  So  does 
the  heat  content  in  the  region  bounding  the  hydrographic  section.  To  which  extent 
estimates  of  heat  content  and  volume  flux  from  one-time  hydrographic  survey  can 
represent  their  mean  values  is  the  issue  to  be  address  here.  As  detailed  later,  the 
method  used  by  Parrilla  et  al.  (1994)  and  Arbic  and  Owens  (1999)  to  estimate  uncer¬ 
tainty  in  average  temperature  (heat  content)  is  empirical.  Here  I  will  provide  a  more 
rigorous  estimate  of  uncertainty  for  the  volume  transport  and  heat  content  (average 
temperature). 

Suppose  that  the  hydrographic  section  runs  in  the  east-west  direction  and  is 
bounded  at  both  ends  with  the  zonal  length  of  L  km.  The  volume  transport  across 
this  section  from  a  depth  of  — /z2  to  ~h\  m  for  each  mode  can  be  written  as 

rL  r  —  h\ 

V(L,y,n,t)  =  /  /  v(x,y,  z,n,t)dzdx ,  (6-2) 

Vo  J— h,2 

and  the  heat  content  in  the  region  bounding  the  hydrographic  section  is 

rL  r—h\ 

Q(L,y,n,t)=  I  /  pocp9(x,y,z,n,t)dzdx.  (6-3) 

Because  the  value  of  the  density  p0  changes  little  with  location  and  time,  here  it  is 
treated  as  a  constant,  po  =  1.027  x  103  kg/m3.  One  important  and  widely  used  factor 
associated  with  heat  content  is  the  average  temperature  in  the  region  bounding  the 
hydrographic  section.  By  definition,  the  average  temperature  can  be  written  as 

1  rL  r-hi 

[e](L,y,n,t)  =  —  /  9(x,  y,  z,  n,t)dzdx,  (6.4) 

tlJ-J  Jo  j  —  /l2 
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where  h  =  /j2  —  h\ .  In  chapter  3,  I  have  derived  an  analytic  formula  for  the  frequency 
spectrum  of  tomographic  data,  equation  (3.37).  From  that  formula,  the  variance  of 
the  average  temperature  can  be  readily  obtained 

9  f  1  f-h  1  (9#n  \  2  Z’+OO  /-H-OO  r+oo 

ofo(L,n,A,0)  =  I -J_h  -^Gn{z)dz J  J  J  J  f2${k,l,u,n,\,<j))W{k,L)dujdkdl. 

(6.5) 

The  variance  of  the  average  meridional  velocity  can  be  obtained  in  a  similar  way: 

/ 1  f~hi  \  2  r+oo  r+oo  r+oo 

a(v](L,n,X,(t>)=  -/  Fn{z)dz  /  /  /  47r2A:2$(A:,/,a;,n,  \,4>)W{k:  L)dcodkdl, 

ytl  J  —h‘2  J  J—  oo  J —  oo  J —o o 

(6.6) 

where  W(k,L )  is  defined  in  equation  (3.38).  Note  that  FF(A;,L)  decreases  with  L 
and  IF(/i:,L  =  0)  =  1.  The  variances  of  the  average  temperature  and  average  merid¬ 
ional  velocity  decrease  with  L  because  more  energy  is  filtered  out  as  L  increases. 
afe]{L,  n,  A,  </>)  and  a fv](L,  n,  A,  </>)  depend  on  location  as  well  as  the  zonal  length  L.  At 
each  site  n ,  A,  4> )  and  cr^(L,  n,  A,  <f>)  are  normalized  by  their  maximum  values, 

°[9]{L  =  0,  n,  A,  4>)  and  cfv](L  =  0,  n,  X,(j>),  so  that  they  only  depend  on  the  zonal 
length  L.  As  displayed  in  figure  6-4,  the  normalized  a^(L,  n,  A,  4>)  and  vfvAL,  n,  A,  <f>) 
decrease  with  the  zonal  length  L.  This  suggests  that  the  longer  the  zonal  length  of  the 
hydrographic  section,  the  less  variable  the  average  temperature  and  mass  transport. 

The  most  striking  property  of  figure  6-4  is  that  a^(L,n,  X,4>)  decreases  with  zonal 
length  L  much  less  rapidly  than  crj^(L,  n,  A,  </>).  For  average  temperature,  half  of  the 
energy  is  filtered  out  when  the  zonal  length  is  about  350  km;  for  average  meridional 
velocity,  90%  of  the  variance  is  removed  due  to  filtering  when  the  zonal  distance  is 
250  km.  This  difference  arises  because  the  zonal  spatial  scale  of  meridional  velocity 
is  shorter  than  that  of  temperature,  as  discussed  in  section  6.1. 

The  variance  of  heat  content  is 


<7%{L,n,\,  <j> )  =  p2<%h2L2ofa(L,n,\,<ti), 


(6.7) 


and  the  variance  of  volume  flux  is 


cry  (L,  n,  A,  <p)  =  h2L2a2v]  (L,  n,X,4>).  (6.8) 

From  equation  (6.7)  and  (6.8),  we  can  see  that  <Jq(L,  n,  A,  0)  and  ay(L,n,  A,(/>)  does 
not  simply  increase  with  L2.  How  <r@(L,  y,  n)  and  &y(L,  y ,  n )  depend  on  L  is  displayed 
in  figure  6-5.  The  variance  of  volume  flux  increases  rapidly  as  L  increases  from  0  to 
200  km,  then  it  increases  very  slowly.  As  L  increases  from  200  km  to  2000  km, 
the  variance  of  volume  flux  only  increases  by  about  10%.  This  phenomenon  can  be 
rationalized  in  the  following  way:  because  the  zonal  length  scale  of  mesoscale  eddies  is 
about  200  km,  the  hydrographic  section  contains  less  than  one  mesoscale  eddy  when 
L  is  shorter  than  200  km,  and  the  variability  will  increase  with  L  until  L  reaches 
200  km.  As  L  is  longer  than  200  km,  the  number  of  mesoscale  eddies  in  the  region 
bounding  the  hydrographic  section  will  be  greater  than  one  and  increases  with  L;  the 
variability  of  total  volume  flux  will  be  the  average  of  variability  of  different  eddies. 
Therefore,  it  changes  slowly  when  the  zonal  length  is  greater  than  200  km.  However, 
because  the  zonal  length  scale  of  temperature  is  much  larger  than  that  of  meridional 
velocity,  the  variance  of  heat  content  still  increases  smoothly  with  L  after  L  is  greater 
than  200  km  (figure  6-5). 

The  variances  of  volume  transport  and  heat  content  depend  on  location  as  well 
as  L.  As  an  example,  the  standard  deviation  of  the  volume  flux  and  heat  content 
for  the  hydrographic  section  at  24°N  is  estimated.  The  zonal  range  of  the  hydro- 
graphic  section  is  assumed  to  be  from  25°W  to  75°W,  which  nearly  coincides  with 
the  actual  hydrographic  section.  The  function  for  l ,  u>,  n,  A,  4>)  along  this  section 

can  be  derived  according  to  the  formulas  produced  in  chapter  5.  For  the  volume 
transport,  I  choose  the  vertical  integration  to  be  from  the  surface  to  the  depth  of 
2000  meters,  while  for  heat  content  (temperature),  the  vertical  integration  is  taken 
from  1000  meters  to  2000  meters  because  the  signal-to-noise  ratio  for  temperature  is 
largest  within  this  range.  The  standard  deviation  of  volume  flux,  average  temperature 
and  heat  content  for  this  particular  hydrographic  section  are  listed  in  table  6-1.  As 
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displayed  in  the  table,  the  variability  of  barotropic  volume  flux  and  baroclinic  volume 
flux  is  roughly  the  same.  For  the  average  temperature,  the  variability  in  the  second 
baroclinic  mode  is  larger  than  that  in  the  first  baroclinic  mode. 


oy  (0-2000  m) 

o[e]  (1000-2000  m) 

a©  (1000-2000  m) 

n=0 

2.9  (Sv) 

n=l 

2.0  (Sv) 

0.01  (°C) 

1.7  x  1014  (J/m) 

n=2 

1.1  (Sv) 

0.018  (°C) 

2.9  x  1014  (J/m) 

total 

3.7  (Sv) 

0.021  (°C) 

3.4  x  1014  (J/m) 

Table  6-1:  Standard  deviation  of  volume  flux,  average  temperature  and  heat  con¬ 
tent  for  the  first  three  modes  and  the  standard  deviation  of  the  total  variability. 

Arbic  and  Owens  (2000)  calculated  the  zonal  average  and  standard  deviation  of 
the  temperature  difference  between  1981  and  1957.  They  assumed  that  the  standard 
error  uncertainty  in  the  zonally  averaged  temperature  difference  is  mainly  due  to 
mesoscale  eddies  that  have  a  horizontal  length  scale  of  about  100  km  and  obtained 
the  error  uncertainty  by  dividing  the  standard  deviation  by  the  square  root  of  the 
number  of  mesoscale  eddy  length  scales  in  the  5,000  km  width  of  the  24°N.  The  zonally 
and  vertically  (from  1000  to  2000  m)  averaged  temperature  difference  between  1984 
and  1957  is  0.11°C  and  the  error  uncertainty  is  0.024°C,  which  is  very  close  to  my 
estimate. 

Ganachaud  (1999)  estimated  the  uncertainty  in  volume  flux  through  an  eddy¬ 
resolving  oceanic  general  circulation  model  and  he  found  that  the  variability  in  baro¬ 
clinic  volume  flux  through  24°  above  2000  m  is  about  4  Sv  (personal  communication), 
which  is  larger  than  the  value  of  2.4  Sv  I  obtained. 
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Zonal  distance  (km) 

Figure  6-5:  The  normalized  <Jq(L,  n,  A,  <f>)  (solid  line)  and  Oy(L,n,  X,4>)  (dashdot 
line)  as  a  function  of  L.  <7Q(L,n,  A,  </>)  and  cry(L,n,  \,(f>)  are  normalized  by  cr©(L  = 
2000,  n,  A,  4>)  and  Oy  (L  =  2000,  n,  A,  <ft),  respectively. 
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6.3  Design  of  an  observational  network 


Fluctuations  of  heat  storage  in  the  ocean  are  key  factors  in  understanding  the  general 
circulation  of  the  ocean  and  climate  change.  Because  the  area  of  the  ocean  is  vast, 
it  is  very  expensive  to  monitor  the  climatic  variability  of  heat  content.  The  oceanic 
variability  is  continuous  in  space  and  time.  The  measurements,  however,  are  sparse 
and  subject  to  noise.  There  is  no  doubt  that  the  higher  the  sampling  density,  the  more 
accurate  the  estimates.  Nonetheless,  there  is  an  inevitable  trade-off  between  cost  and 
sampling  density.  It  is  impossible  and  unnecessary  to  make  the  sampling  spacing  in¬ 
finitesimal.  It  is  desirable  to  evaluate  the  effectiveness  of  different  sampling  strategies 
before  the  monitoring  network  is  put  into  use.  The  particular  question  posed  here 
is:  to  what  extent  can  the  finite  and  discrete  observations  can  represent  the  actual 
oceanic  variability?  One  closely  related  question  is:  what  is  the  minimum  sampling 
density  in  order  to  resolve  and  detect  the  climate  change  with  certain  accuracy? 

6.3.1  Optimal  interpolation 

Optimal  estimation  methods  have  been  applied  by  a  few  authors  to  the  design  of 
observational  networks  in  oceanography  (e.g.,  Bretherton  et  ah,  1976;  White  and 
Bernstein,  1979;  Bretherton  et  al,  1984;  White,  1995).  The  advantage  of  optimal 
estimation  theory  is  that  one  can  quantitatively  estimate  the  impact  of  sampling 
strategy  and  measurement  noise  on  the  final  product.  An  important  factor  in  ap¬ 
plying  this  theory  is  the  knowledge  of  the  statistical  structure  of  the  variable  to  be 
analyzed.  Once  the  autocorrelation  and  signal-to-noise  ratio  are  determined,  then 
the  least  sampling  rates  that  produce  the  desirable  interpolation  error  can  be  esti¬ 
mated.  The  general  approach  to  estimate  the  random  vector  X  from  its  measurements 
Y  in  statistical  estimation  theory  is  the  so-called  Bayesian  least-squares  estimation 
method  (Papoulis,  1991).  The  key  factor  of  this  approach  is  the  conditional  proba¬ 
bility  px|y(X|Y),  which  fully  characterizes  the  relationship  between  the  signal  vector 
X  and  the  observation  vector  Y.  In  general,  this  estimator  is  a  nonlinear  function 
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of  the  data  and  it  yields  the  least  uncertainty  (in  terms  of  error  variance)  among  all 
estimators  (Willsky  et  ah,  1998).  The  difficulty  inherent  in  this  method  is  the  deriva¬ 
tion  of  the  conditional  probability  px|y(X|Y),  which  often  is  not  available.  Another 
widely  used  approach  is  minimum  variance  (Gauss-Markov)  estimator.  By  assuming 
that  the  estimator  to  be  a  linear  function  of  the  data,  this  estimator  requires  only 
the  first  and  second  statistical  moments:  mean  and  covariance  distribution  of  the 
signal  and  noise.  The  minimum  variance  estimation  method  was  first  introduced  into 
meteorology  by  Gandin  (1963)  to  estimate  the  value  of  the  field  at  a  point  based  on 
randomly  distributed  measurements.  Its  application  to  oceanographic  problems  was 
reviewed  by  Wunsch  (1996).  Assuming  that  the  field  and  noise  are  jointly  Gaussian 
distributed,  Bretherton  et  al.  (1984)  applied  the  Bayesian  least-squares  estimation 
method  to  the  design  study  of  measurements  of  heat  storage.  Under  this  assump¬ 
tion,  the  Bayesian  least-squares  method  and  linear  least-squares  method  produce  the 
same  estimate  and  the  same  uncertainty  (Willsky  et  al.,  1998).  In  the  following  I 
will  use  the  linear  minimum  variance  method  to  evaluate  the  effectiveness  of  different 
sampling  strategies. 

The  minimum  variance  estimator 

Let  rj  =  ri,  r2, ...,  tk  denote  the  locations  where  the  measurements  are  made,  and 
let  sm  =  Si,  s2, ...,  sl  denote  the  positions  where  the  estimates  are  sought  and  which 
are  usually  at  grid  point.  Define  the  measurement  at  the  point  ri  as  y(ri),  which  is 
the  actual  value  plus  some  random  noise 

y(r,)  =  x(ri)  -f  w( r;),  1  <  i  <  K.  (6.9) 

I  assume  that  the  noise  has  zero  mean,  the  noise  at  different  points  is  uncorrelated 
and  the  noise  is  uncorrelated  with  the  actual  value  of  the  underlying  process;  that  is 

<  w(t{)  >=  0,  <  w(ri)w(rj)  >=  <4(ri)<%,  <  (x(rj)  -  x(ri ))w(rj)  >-  0,  (6.10) 

where  x(ri)  represents  the  mean  value  of  x(r5). 
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The  estimate  at  the  location  sm  is  assumed  to  be  an  linear  combination  of  the 
measurements, 

K 

^(Sm)  G'(Sm)  T  ^  ri)j/(ri)-  (6.11) 

i=  1 

In  order  for  the  estimate  to  be  unbiased,  I  choose 

K 

a(sm)  =  ^(r;)  -  Y,  K sm,  ri)y(rj).  (6.12) 

i=  1 

The  weighting  coefficient  6(sm,rj),  which  minimizes  the  mean-square  error,  is  deter¬ 
mined  by  the  following  simultaneous  equations 

K 

^(®m?  rj)pix(rij  rj)  +  7  KSm?  ri)  =  Pxx(sm7  ri)j  1  ^  ^  (6.13) 

3= 1 

where  Pzx^rj)  is  the  autocorrelation  between  signals  at  the  location  ri  and  rj, 
Pxx  (®m?  Ti)  is  the  signal  autocorrelation  between  the  point  where  the  interpolation 
takes  place  and  the  surrounding  points  where  observations  are  made,  and  y2  is  the 
ratio  of  signal  variance  to  noise  variance. 


.,2  _  °x 


7  = 


at 


The  corresponding  mean-square  error  is 


K 


2=1 


(®m)  CTX[!  ^i)Pxx  (Sm?  ^l)] . 

Define  the  normalized  error  as 


(6.14) 


(6.15) 


therefore 


® e  (^m) 

5 

°x 


P  = 


1 


K 


2=1 


(6.16) 


(6.17) 
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Equation  (6.17)  was  derived  by  Alaka  and  Elvander  (1972)  and  used  by  White  and 
Bernstein  (1979)  to  design  an  oceanographic  network.  Equation  (6.17)  states  that 
the  interpolation  error  depends  only  on  the  autocorrelation  function,  interpolation 
scheme  and  weighting  coefficients.  It  is  independent  of  the  values  of  measurements. 
Therefore,  the  interpolation  error  can  be  estimated  before  any  measurements  are 
made  as  long  as  the  autocorrelation  function  and  error  statistics  are  known. 

The  correlation  functions  for  temperature  and  two  components  of  horizontal  ve¬ 
locity  have  been  obtained  in  section  6.1.  The  signal  is  often  defined  by  the  scales  of 
interest  to  the  observer  and  usually  stands  for  large  scale  perturbations.  The  noise, 
on  the  other  hand,  has  two  sources:  one  is  the  instrumental  error,  and  the  other  is 
the  subgrid  ambient  variability  in  the  field  itself.  In  the  literature,  the  noise  is  often 
considered  to  be  dominated  by  the  second  source.  Once  the  wavenumber  spectrum 
of  field  and  sampling  spacing  are  known,  the  variance  of  noise  due  to  subgrid  aliasing 
can  be  easily  calculated.  Based  on  the  zonal-wavenumber  spectrum  of  temperature 
obtained  in  last  chapter,  it  is  found  that  41%  of  the  variance  is  distributed  at  zonal- 
wavelengths  longer  than  1000  km,  64%  at  zonal-wavelengths  longer  than  500  km  and 
92%  at  zonal-wavelengths  longer  than  200  km.  Therefore,  the  ratio  of  signal  vari¬ 
ance  to  error  variance  is  about  0.7:1  for  the  sampling  array  with  the  zonal  spacing 
of  500  km.  It  increases  to  1.8:1  for  zonal  spacing  of  250  km  and  to  12:1  for  the 
zonal  spacing  of  100  km.  The  above  discussion  is  about  temperature.  The  results 
for  horizontal  velocities  are  different.  Because  the  zonal-wavenumber  spectrum  of 
temperature  is  redder  than  that  of  meridional  component  of  horizontal  velocities,  the 
ratio  of  zonal-wavenumber  spectrum  of  temperature  to  the  zonal-wavenumber  spec¬ 
trum  of  meridional  velocity  is  proportional  to  k~ 2.  It  is  found  that  only  10%  of  the 
variance  of  meridional  velocity  is  at  zonal-wavelengths  longer  than  500  km  and  33  % 
at  wavelengths  longer  than  200  km.  Therefore,  for  the  meridional  velocity,  the  ratio 
of  signal  variance  to  error  variance  is  about  0.11:1  for  the  sampling  array  with  the 
zonal  spacing  of  250  km  and  0.5:1  for  the  zonal  spacing  of  100  km.  The  ratio  of  signal 
variance  to  error  variance  for  temperature  is  much  higher  than  that  for  meridional 
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velocity  for  the  sampling  array  with  zonal  spacing  larger  than  100  km.  Here  I  am  not 
going  to  specify  the  value  of  the  signal-to-noise  ratio.  Instead,  I  let  it  to  be  a  control¬ 
lable  variable,  and  I  will  study  how  the  interpolation  error  depends  on  signal-to-noise 
ratio  as  well  as  sampling  density. 

In  the  most  general  cases,  the  interpolation  is  four-dimensional  in  zonal,  merid¬ 
ional,  temporal  and  vertical  directions.  Here,  I  consider  the  interpolation  to  be  two- 
dimensional  and  in  the  horizontal  plane.  In  this  case,  the  autocorrelation  function 
is  two-dimensional:  Rg(rx,ry )  =  Rg{rx)Rg(ry),  where  Re{rx )  and  Rg{ry)  are  the  one¬ 
dimensional  zonal  autocorrelation  function  and  meridional  autocorrelation  function, 
respectively,  and  have  been  obtained  in  section  6.1.  White  (1995)  compared  the 
results  in  different  dimensions  and  found  that  the  interpolation  error  decreases  as 
dimension  increases  because  there  are  greater  independent  observations  for  higher 
dimensions. 

As  an  illustration  of  how  the  normalized  interpolation  error  varies  with  sampling 
density,  I  choose  a  12-point  interpolation  scheme,  which  is  displayed  in  figure  6-6. 
I  want  to  estimate  the  temperature  at  the  center  point  based  on  measurements  at 
the  surrounding  12  points.  As  displayed  in  figure  6-7,  the  interpolation  error  drops 
as  sampling  density  increases,  more  rapidly  in  the  beginning.  The  less  accurate  the 
observations,  the  more  the  interpolation  error.  The  structure  of  mapping  errors  in 
figure  6-7  suggests  that  optimal  sampling  occurs  at  spacing  about  40-50  km,  because 
interpolation  error  decreases  rapidly  for  sampling  spacing  greater  than  50  km,  then 
very  slowly. 

As  mentioned  before,  my  model  is  a  first-order  approximation  to  the  actual  ocean 
variability.  The  dependence  of  spectral  shape  on  geography  has  been  simply  neglected. 
However,  as  described  by  Stammer  (1997),  the  observed  spectral  shape  does  depend 
on  geography  and  latitude.  Therefore,  for  a  more  precise  design  for  sampling  some 
particular  regions,  the  particular  statistics  of  those  regions  should  be  used.  Figure  6-8 
shows  how  sensitive  the  interpolation  error  is  to  the  spectral  shape.  The  three  curves 
in  figure  6-8  display  a  similar  structure,  which  suggests  that  the  optimal  sampling 
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spacing  for  the  three  different  spectral  shapes  is  roughly  the  same,  around  40-50  km. 
When  the  spectral  slope  changes  by  15%,  the  variations  of  the  normalized  error  is 
less  than  0.1. 

The  strategy  for  monitoring  horizontal  velocities  should  be  different  from  the 
strategy  for  measuring  temperature,  because  horizontal  velocities  have  distinct  auto¬ 
covariance  structure.  The  most  important  property  of  the  auto-covariance  of  hor¬ 
izontal  velocities  is  anisotropy;  the  zonal  decorrelation  scales  are  longer  than  the 
meridional  scales  for  zonal  component,  and  the  reverse  is  true  for  meridional  compo¬ 
nent. 

There  is  an  emerging  program,  named  Argo,  which  is  going  to  deploy  a  global 
array  of  3,000  profiling  floats  to  measure  the  temperature  and  salinity  of  the  ocean’s 
upper  layer.  If  the  3,000  floats  are  deployed  regularly  in  the  global  ocean,  the  spacing 
will  be  about  340  km  by  340  km.  The  wavenumber  spectrum  of  temperature  suggests 
that  56%  of  the  variance  is  at  wavelengths  longer  than  680  km;  therefore,  the  ratio 
of  signal  variance  to  noise  variance,  y2,  is  about  1.3  for  the  sampling  array  with  340 
km  spacing.  Figure  6-7  shows  that  when  the  ratio  of  signal  variance  to  noise  variance 
is  1.5,  50  km  sampling  spacing  is  preferable.  If  we  further  double  the  sampling 
density — sampling  spacing  is  decreased  from  50  km  to  25  km — the  interpolation  error 
will  decrease  only  10%.  If  the  floats  are  regularly  distributed  with  spacing  50  km, 
about  140,000  floats  are  needed  to  cover  the  global  ocean.  Therefore,  the  proposed 
3000  floats  cannot  monitor  the  dominant  variability  over  the  global  ocean. 

6.3.2  Optimal  averaging 

In  the  above,  I  studied  how  to  estimate  the  value  at  a  point,  based  on  some  observa¬ 
tions  surrounding  that  point.  For  the  study  of  climate  change,  it  is  highly  desirable 
to  obtain  estimates  of  the  value  averaged  over  a  large  area,  rather  than  the  value  at 
a  point.  Kagan  (1997)  gave  a  comprehensive  and  detailed  discussion  of  statistical 
averaging  of  meteorological  fields. 


99 


o 


Figure  6-6:  The  interpolation  scheme  used  in  this  study.  The  solid  circles  represent 
the  places  where  the  measurements  are  made  and  the  open  circle  stands  for  the 
location  where  the  estimate  to  be  obtained. 


station  spacing  (km) 


Figure  6-7:  Normalized  interpolation  error  of  temperature  versus  sampling  density 
for  different  values  of  y2,  the  ratio  of  signal  variance  to  noise  variance. 
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Figure  6-8:  The  normalized  error  for  different  wavenumber  spectral  shapes  with 
signal-to-noise  ratio  of  1.0.  The  solid  line  corresponds  to  the  globally  averaged 
wavenumber  spectral  shape,  while  the  dashdot  and  dotted  line  correspond  to  the 
wavenumber  spectrum  with  slope  15%  more  and  15%  less  than  that  of  the  globally 
averaged  one,  respectively. 

In  the  following,  the  optimal  averaging  method  will  be  applied  to  the  design  of 
oceanographic  monitoring  network,  with  a  focus  on  studying  how  the  error  associated 
with  determining  averaged  value  depends  on  sampling  density.  I  assume  that  the 
averaging  is  two-dimensional  and  in  the  horizontal  plane.  If  F(x,  y)  is  defined  to  be 
the  averaged  value  of  the  field  f(x,y )  (figure  6-9),  then 


1  rxs+Lx/2  rys+Ly/2 

F(x„y,)  =  1-r  /  /(fi,6)de 

-Ls 'Y' L/ 7#  «/ jPo  L t  / 2  J xi a  —Li,  12 


Jx*-Jy  ^Xs  Lx/ 2  Jys  ^y/2 


(6.18) 


Let  g(xi,  yi)  be  the  measurement  at  the  point  ( Xi,yi ),  which  is  the  true  value  f(x{,yi) 
plus  some  noise  w(xi,yi),  so  that 


g{xi,  yi)  =  f(xi,  yi)  +  wfa,  yi ). 


(6.19) 
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The  estimate  of  F(x,  y )  is  assumed  to  be  an  linear  combination  of  surrounding  mea¬ 
surements: 

n 

F{xs,ys)  =  Y.P*9(xi,Vi),  (6.20) 

i= 1 

where  pi  is  the  weighting  factor  for  the  measurement  at  the  point  ( Xi ,  y,). 

For  a  homogeneous  process,  in  order  to  get  unbiased  estimate,  pi  needs  to  satisfy 

n 

YlPi  = (6-21) 

2= 1 

unless  <  /  >=  0.  In  this  study,  I  assume  that  the  time-mean  has  been  removed 
before  the  spatial  average  is  made.  I  further  assume  that  actual  value  is  uncorrelated 
with  noise,  noises  at  different  places  are  uncorrelated  and  that  noises  have  zero  mean. 
The  mean-square  error  is 


E2{xs,  ys)  =<  [F(xs,  ys)  -  F(xs,  ys )]2  >  .  (6.22) 


For  an  unbiased  estimate, 

E  (:rs,  ys )  —  —  2  ^2,PiRf f.  +  PiPjFjj  +  &w  ^2  Pi>  (6.23) 

2=1  2=1  J  =  1  2=  1 

where  o\  is  the  variance  of  the  average  value,  Rtj  is  the  covariance  between  signals 
at  the  point  ( a^y* )  and  ( Xj,yj ),  a £  is  variance  of  noise  and  RFfi  is  the  covariance 
between  the  average  value  and  the  value  at  the  point  ( Xi ,  y,), 


RFfi  -<  [f(xi,yi)~  <  >][F(xs,ys)-  <  F(xs,ys )  >]  >  .  (6.24) 


Substitution  of  equation  (6.18)  into  (6.24)  yields 


1  rxs-\-Lx  {2  rys-\-Ly  /2 

rfu  =  yt~  /  ,  /  Rf(x*  ~  &’Vi  ~  &)dtidt2, 

J x$  Lx/2  J ys — Z/y/2 


(6.25) 


where  Rf(xi  —  ,  y*  —  £2)  is  covariance  between  values  at  point  ( o^y, )  and  (£l5£2). 
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Equation  (6.25)  suggests  that  covariance  between  the  average  value  and  the  value  at 
the  point  ( Xi,yi )  is  the  average  of  the  covariances  between  the  value  at  point  {xi^yi) 
and  at  every  point  inside  the  area  where  the  average  to  be  estimated. 

For  simplicity,  define  (rx,  ry )  to  be  the  separation  distance  between  the  point 
and  the  center  of  the  area  where  the  average  value  to  be  estimated, 


rx  =  x i  -  xs,  ry  =  y{  —  ys, 


(6.26) 


so  that  equation  (6.25)  can  be  simplified  as 


Rfu  = 


1  [Tx  +  Lx/2  rry+Ly/2 


/Tx+bx//  r 

J  v x.  Lx  / 2  J v i 


^x-^y  tx~Lx f 2  Lyf2 


Rf(Tl,T2)dT1dT2. 


(6.27) 


Normalize  (6.23)  by  the  variance  of  field  f,cr2, 


2  E2(xs,ys)  ^  0  ,  ^  ,  _2^  2 

e  = - o - =\i-  2}_^PiPi  +  2^2^PiPjrij  +  7  2^Pzi 


o. 


(6.28) 


2—1 


2=1  j=  1 


2=1 


where 


/*  = 


A  = 


R 


Ffi 


°2f 


(6.29) 


Note  that  the  variance  of  the  average  value  is  not  equal  to  the  average  of  the  variance 
at  single  point.  In  fact  the  variance  of  the  average  value  is  always  less  than  variance  of 
the  value  at  a  point,  p  <  1,  because  variability  associated  with  small  scales  is  filtered 
out  due  to  averaging.  For  temperature,  the  variance  at  a  point  is 


/-f-OO  /’-f-00  r-\~  00 

/  /  \@a\2&(k,l, oj,  n,  \,<j))dkdldu, 

-OO  7—00  7  —  00 


(6.30) 


and  the  variance  of  the  average  value  is 


/-Too  i — boo  /•'Too 
oo  7—oo  7—oo 


sin2(2irkLx/2)  sin2(2'KlLy/2 ) 
{2MJW  (2nlLy/2)2 


dkdldu 

(6.31) 
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If  Lx  =  1000  km  and  Ly  =  1000  km,  then  p  =  0.05. 

The  optimal  weight  is  obtained  by  minimizing  the  error  variance  e2.  Therefore,  p 
can  be  obtained  by  solving 

n 

Pjrij  +  Pil~2  =  A,  (6.32) 

j~  1 

and  the  corresponding  error  is 


n 

=  Pifc-  (6.33) 

2=1 

The  error  depends  on  the  size  of  the  area  where  the  average  value  to  be  estimated,  the 
number  of  observations,  the  sampling  density  and  signal-to-noise  ratio.  The  following 
discussions  describe  two  simple  cases.  First,  the  averaging  area  is  fixed  to  be  1000 
km  by  1000  km,  the  number  of  observations  is  fixed  to  be  12,  and  the  average  scheme 
is  shown  in  figure  6-9.  Figure  6-10  shows  how  the  error  depends  on  sampling  spacing 
and  signal-to-noise  ratio.  As  shown  there,  the  normalized  error  is  always  less  than 
0.22,  which  is  the  square  root  of  the  normalized  variance  of  the  average  value.  This 
means  that  variance  of  the  error  cannot  exceed  the  variance  of  the  average  value.  The 
most  striking  property  in  figure  6-8  is  that  the  minimum  errors  are  obtained  at  the 
sampling  spacing  around  220  km.  This  is  because  the  number  of  observations  is  fixed 
to  be  12.  When  the  measurements  are  too  close  to  each  other,  they  just  duplicate  one 
another.  Some  improvements  can  be  achieved  by  mutual  canceling  of  random  error 
in  the  data  by  repeated  sampling.  If  they  are  too  far  away  from  each  other,  however, 
they  will  miss  the  dominant  variability  of  the  underlying  process.  In  the  second  case, 
the  averaging  area  is  also  fixed  to  be  1000  km  by  1000  km.  However,  the  number  of 
observations  is  not  fixed,  and  I  assume  that  all  observations  are  located  within  the  area 
where  the  average  value  is  estimated  (figure  11).  As  shown  in  figure  6-12,  the  error 
associated  with  determining  the  average  value  decays  smoothly  with  sampling  density 
and  more  rapidly  at  sampling  density  less  than  0.012  station/km;  and  it  is  much  less 
than  the  error  associated  with  estimating  the  value  at  a  point  (figure  6-7).  Because 
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the  root  mean  square  of  the  average  value  is  about  0.22,  the  error  associated  with 
determining  the  average  value  is  always  less  than  0.22.  Both  of  the  two  cases  studied 
here  have  limitations.  In  the  first  one,  the  number  of  observations  used  to  estimate 
the  average  value  is  fixed  to  12;  in  the  second  case  no  observations  outside  the  average 
area  have  been  used  to  estimate  the  average  value.  In  order  to  fully  understand  the 
relationship  among  the  error,  samping  density,  number  of  observations,  and  sampling 
geometry,  further  study  is  needed. 


6.4  Lagrangian  correlation 

Recently,  there  have  been  an  increasing  number  of  observations  from  drifting  floats. 
It  is  of  interest  to  find  out  how  to  combine  float  data  with  other  measurements  such  as 
XBT  and  altimeter  data  for  the  purpose  of  understanding  the  ocean  general  circula¬ 
tion.  The  relationship  between  Lagrangian  and  Eulerian  statistics  is  fundamental  to 
the  understanding  of  turbulent  diffusion.  A  number  of  theories  claim  to  give  the  re¬ 
lationship  between  the  velocity  statistics  of  Lagrangian  particles  and  the  statistics  of 
the  Eulerian  flow  (see  the  review  by  Davis,  1981).  These  theories  enable  one  to  relate 
the  Lagrangian  autocovariance  function  to  the  Eulerian  energy  wavenumber  and  fre¬ 
quency  spectra.  Therefore,  given  Eulerian  energy  wavenumber  and  frequency  spectra, 
one  can  predict  the  Lagrangian  autocovariance  function  (frequency  spectrum)  based 
on  various  theories.  The  crucial  factor  to  test  these  theories  is  the  three-dimensional 
spectrum  of  the  Eulerian  flow.  The  motivation  here  is  to  test  one  of  the  theories  by 
comparing  the  predicted  Lagrangian  autocorrelation  function  with  that  of  observa¬ 
tions.  One  widely  used  approximation  for  the  relationship  between  Lagrangian  and 
Eulerian  statistics  in  the  literature  was  originally  proposed  by  Corrsin  (1959)  and 
applied  by  Saffman  (1963).  The  theory  by  Corrsin  will  be  evaluated  here. 

I  assume  that  the  particle  trajectory  is  in  the  horizontal  plane  and  let  (X’(i),  Y(t)) 
be  the  trajectory  of  the  particle,  which  was  deployed  at  the  position  of  (X(t)  =  0, 
Y ( t )  =  0)  at  time  t  =  0;  then  the  trajectory  satisfies  the  equations 
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Figure  6-9:  The  optimal  averaging  scheme  used  in  this  study.  The  solid  circles 
represent  the  places  where  the  measurements  are  made  and  the  open  circle  stands  for 
the  center  of  the  area  where  the  average  are  performed. 


Figure  6-10:  Normalized  interpolation  error  of  temperature  versus  sampling  spac¬ 
ing  for  different  values  of  72,  the  ratio  of  signal  variance  to  noise  variance. 
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Figure  6-11:  The  optimal  averaging  scheme  for  the  second  case.  Observations  are 
all  located  inside  the  area  where  the  average  is  made  and  are  regularly  distributed. 
The  spacing  will  decrease  as  the  number  of  observations  increases. 


sampling  density  (station/km) 

Figure  6-12:  Normalized  interpolation  error  of  temperature  versus  sampling  spac¬ 
ing  for  different  values  of  72,  the  ratio  of  signal  variance  to  noise  variance. 


(6.34) 

^M=VL(t)=v(X(t),Y(t),t), 

(6.35) 

where  ui(t )  and  VlO)  are  zonal  and  meridional  components  of  Lagrangian  horizontal 
velocities,  and  where  u(X(t),Y(t),t)  and  v(X(t),Y(t),t)  are  zonal  and  meridional 

Eulerian  horizontal  velocities  at  the  position  ( X(t ),  T(t))  at  time  t.  The  autocovari¬ 
ance  function  of  the  Lagrangian  velocity  ui(t)  is  defined  as 

RLu{t)  =<  uL(0)uL{t)  >=<  u(0,0,0)u(X(t),Y(t),t)  >  . 

(6.36) 

The  particle  displacement  in  the  direction  of  x-axis  at  time  t  is 

X(t)  =  f  uAt)dt . 

Jo 

(6.37) 

Let  a2(t )  be  the  mean  square  particle  displacement  in  the  direction  of  x-axis,  then 

a2(t)  =<  X2{t)  >  . 

(6.38) 

It  is  readily  proved  that 

a2(t)  =  2  [  (t  -  t)RLu(t)(1t. 

Jo 

(6.39) 

Taking  the  second  derivative  of  equation  (6.39)  gives 

d2o2(t)  nn  /x 
d2t  ~  2Rlu^' 

(6.40) 

Equation  (6.36)  can  be  written  as 

r+oo  r+ oo 

RluO)  =  /  <  u( 0, 0,  0)u(rx,  ry ,  t )  S  [rx  -  X (t),  ry  —  Y (t)]  >  drxdry, 

J  —  OO  J  —  oo 

(6.41) 
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where  5[rx  —  X(t),ry  —  Y(t) ]  is  a  two-dimensional  delta  function.  The  reader  is 
reminded  that  [x(t),y(t)]  is  the  position  of  the  particle  at  time  t  and  is  random, 
while  [rx,  ry]  is  the  position  of  a  point  in  two-dimensional  space  and  is  determin¬ 
istic.  5  [rx  -  X(t),ry  -  Y(t)]  =  0  except  when  X(t)  —  rx  and  Y(t)  =  ry.  So 
5[rx  -  X(t),ry  -Y(t)]  =  0  is  the  marker  of  a  fluid  particle.  For  convenience,  de¬ 
fine  a  new  variable  as 

P(r„  r„  <)  =<  6  [r,  -  X(t),  r,  -  Y(t)\  >,  (6.42) 

which  stands  for  the  probability  density  of  the  position  of  a  particle  at  time  t. 

If  we  assume  that  w(0, 0, 0)u(rx,ry,t)  and  6  [rx  —  X(t),ry  —  Y(t),t]  in  equation 
(6.41)  are  statistically  independent,  equation  (6.41)  becomes 

/+oo  f+ oo 

/  Ru(rx,ry,t)P(rx,ry,t)drxdry.  (6.43) 

-OO  J  —  OO 

This  is  the  independence  approximation  suggested  by  Corrsin  (1959).  The  essence 
of  the  independence  approximation  is  that  after  a  long  time,  the  diffusion  process 
is  effectively  the  same  as  a  random  walk.  For  a  random  walk,  the  local  velocity 
and  the  displacement  of  a  particle  after  a  few  steps  are  uncorrelated.  Weinstock 
(1976)  investigated  the  validity  of  this  approximation.  He  showed  that  it  is  the  first 
term  in  a  systematic  expansion  of  the  Lagrangian  function  and  that  the  independence 
approximation  is  satisfied  for  homogeneous  turbulence.  Equation  (6.43)  suggests  that 
the  Lagrangian  autocovariance  is  the  average  of  the  Eulerian  autocovariance,  and  the 
average  is  taken  with  respect  to  the  probability  of  the  position  of  the  particle. 

The  distribution  of  P(rx,ry,t)  is  assumed  to  be  isotropic  and  normal 


Substitution  of  equation  (6.44)  into  (6.43)  and  use  of  Parseval’s  theorem  yield 


/-(-Oo  /* -f-oo 

/  ru(k,l,t)e~2^k2+l2^dkdl, 

-oo  J  — oo 


(6.45) 


where 


n  +  OO 

Ru(rx,ry,t)e2^kr^drxdry 

-OO 


(6.46) 


Because  X(0)  =  0  and  a2(t)  =<  X2(t )  >,  the  initial  conditions  are 


a2(t)  =  0, 


da2(t) 


=  0  at  t  =  0. 


(6.47) 


The  reader  is  reminded  that  the  Eulerian  covariance  function  is 


/+oo  /*+ oo 

/  T  u(kJ,t)dkdL 

-oo  J  —  oo 


(6.48) 


It  follows  from  equation  (6.45)  that 


\RLu(t)\ 


=  |-Ru(t)|  if  t  =  0 
<  \Ru(t)\  if  t  ^  0  , 


(6.49) 


which  shows  that  the  Lagrangian  autocorrelation  is  less  than  the  Eulerian  autocorre¬ 
lation  except  at  t  =  0. 


Substitution  of  equation  (6.45)  into  (6.40)  yields 


d2a2{t) 


/4- oo  /’-}-00 

/  T  u(k,l,t)e-2n2^2+l2^dkdl. 

-oo  J  —  oo 


(6.50) 


This  is  an  integro-differential  equation  for  o2{t)  in  terms  of  the  Eulerian  spectrum 
function  T „(£;,  /,  t).  The  formula  for  T u(k,  l,  t)  can  be  derived  from  chapter  5, 

z,t)  =  E0(n)I((j),  \)F2(z)Bn(k)47r2l2Cn(l)  j  Dn(u>)cos(2nujt)du  (6.51) 

J  —  OO 
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Given  T u(k,  l,  t),  equation  (6.50)  can  be  solved  numerically  with  the  initial  conditions 
given  in  equation  (6.47).  Lundgren  and  Pointin  (1976)  assumed  some  idealized  forms 
for  Yu(k,l,t )  and  estimated  the  Lagrangian  velocity  autocorrelation  function.  They 
found  that  the  results  based  on  the  independence  approximation  by  Corrsin  (1959) 
are  consistent  with  numerical  simulations.  Because  equation  (6.50)  is  nonlinear,  the 
solution  depends  not  only  on  the  spectral  shape  but  also  the  energy  level.  As  an 
example,  I  choose  the  energy  magnitude  to  be  120  cm2/s2,  which  is  a  typical  value 
in  the  mid-ocean.  The  Lagrangian  and  Eulerian  autocorrelation  functions  of  the  first 
baroclinic  mode  at  the  sea  surface  are  displayed  in  figure  6-13.  As  displayed  there, 
both  Eulerian  and  Lagrangian  autocorrelation  functions  decay  exponentially  with 
time  lag.  The  most  striking  feature  in  figure  6-13  is  that  the  Lagrangian  autocor¬ 
relation  is  much  smaller  than  the  Eulerian  autocorrelation.  The  e-folding  temporal 
scales  of  Lagrangian  and  Eulerian  correlation  are  4  and  40  days,  respectively.  There 
are  very  few  observational  studies  in  the  literature  about  the  relationship  between 
Eulerian  and  Lagrangian  statistics.  Using  the  Lagrangian  correlation  predicator  from 
Corrsin  (1959),  Middleton  (1985)  explored  in  detail  the  relationships  between  La¬ 
grangian  statistics  and  Eulerian  frequency  and  wavenumber  spectra.  He  found  that 
in  general  Lagrangian  statistics  are  insensitive  to  the  details  of  Eulerian  frequency 
and  wavenumber  spectrum.  However,  Lagrangian  frequency  spectrum  depends  on  the 
integral  length  and  time  scale  of  the  Eulerian  velocity  field.  Davis  (1985)  analyzed  the 
drifter  measurements  made  off  the  Northern  California  coast  and  found  that  the  La¬ 
grangian  time  scale  is  significantly  shorter  than  the  Eulerian  time  scale.  This  suggests 
that  the  independence  approximation  theory  is  consistent  with  the  observations. 


Ill 


Chapter  7 


Energy  distribution  in  k,  l  and  u 
space  part  II:  a  directional  form 

7.1  Fitting  l,  lj,  n,  </>,  A)  from  observations 

One  of  the  important  properties  of  low  frequency  variability  is  anisotropy:  the  energy 
of  the  sea  surface  height  is  not  evenly  distributed  traveling  in  all  directions.  As  shown 
in  chapter  4,  the  energy  of  westward-going  motions  is  higher  than  that  of  eastward¬ 
going  motions  by  a  factor  of  three  at  low  frequencies  and  large  scales  in  the  area  of 
25°-50°N  and  195°-225°E.  However,  as  will  be  seen  in  chapter  8,  the  partition  of  energy 
among  eastward-going,  westward-going  and  standing  motions  depends  on  geography. 
In  general,  the  energy  associated  with  westward-going  motions  decreases  with  latitude 
and  there  is  more  energy  going  westward  west  of  major  topographic  features.  In  this 
chapter,  I  am  trying  to  model  the  directional  property  of  the  observed  spectrum 
described.  This  model  is  not  universal:  it  only  applies  to  the  particular  region  in 
the  North  Pacific  ocean  where  the  full  three-dimensional  spectra  were  analyzed.  As 
shown  in  chapter  4,  the  difference  between  the  energy  going  southward  and  northward 
is  so  little  that  it  is  simply  ignored  here.  Because  the  frequency  and  wavenumber 
separate  form,  equation  (4.3),  cannot  provide  any  information  about  the  directional 
property  of  the  spectrum,  a  frequency/zonal-wavenumber  coupled  form  is  chosen  here 
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to  differentiate  the  energy  of  eastward-going  motions  from  that  of  westward-going 
motions.  Again,  because  the  spectrum  is  symmetrical,  in  the  sense  that  T(k,l,cu)  — 
am  only  concerned  with  positive  frequencies,  while  k  and  l  have 
both  signs.  Equation  (4.3)  is  now  modified  to 


$(k,l,u,n,  cf>,  X)  =  Un{k,u)Cn(l)Eo(n)I((j>,  X). 


(7.1) 


nn(ft,cu)  has  distinct  forms  for  eastward-going,  westward-going  and  standing  motions, 


n„(fc,  w) 


Qn—wBn—w(k^Dn_w  (cu) 
<  On-e^n-e(^)^n-e(w) 

Qn— O^n—O  (k)Dn_0  (tu) 


for  k  <  —k0,u  >  0 
for  k  >  ko,uj  >  0 
for  -  k0  <  k  <  k0,  u  >  0 


(7.2) 


where  Bn_w(k),  Bn_e(k )  and  Bn_0(k)  represent  the  wavenumber  spectral  shape  of 
eastward-going,  westward-going  and  standing  motions,  respectively;  Dn_w{uj),  Dn_e(uj) 
and  Dn_0(uj)  stand  for  the  corresponding  frequency  spectral  shape;  and  Qn-W,  Qn-e 
and  Qn- o  determine  the  relative  energy  level  of  eastward-going,  westward-going  and 
standing  motions.  Theoretically,  the  wavenumber  of  standing  motions  is  zero  and 
the  energy  distribution  for  standing  motions  is  a  delta  function  at  k  =  0.  However, 
the  spatial  length  of  the  actual  data  is  finite;  therefore,  any  motion  with  wavelengths 
longer  than  the  spatial  length  of  the  data  can  not  be  resolved  and  is  treated  as  stand¬ 
ing  waves.  In  order  to  be  consistent  with  observations,  I  define  &o  as  the  smallest 
wavenumber  which  can  be  resolved  by  the  data  and  treat  the  motions  with  wavenum¬ 
bers  less  than  k0  as  standing  waves  in  my  model.  For  the  case  analyzed  in  chapter  4, 
k0  is  about  0.0004  cycles  per  km. 

Choosing  the  n„(A;,a;)  to  fit  the  observed  zonal-wavenumber/frequency  spectrum 
of  sea  surface  height 

For  westward-going  motions  ( k  <  -k0,u  >  0) 
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Choosing  the  Bn~w(k)  to  fit  the  observed  zonal-wavenumber  spectrum  of  westward¬ 
going  motions 


Bn-w(k)  \ 


ai|fc|  0-5 

\k\~4 


if  —  hi  <  k  <  —ko 
if  k  <  k\  , 


(7.3) 


where  hi  =  0.0028  (cycles  per  km).  For  continuity  of  Bn( k)  at  k\,  aq  =  3'5. 

Equation  (7.3)  gives 


Bn-W(k)dk  =  7.2  x  107. 


(7.4) 


Choosing  the  Dn-W(u>)  to  fit  the  observed  frequency  spectrum  of  westward- going 
motions 


Bn—w{p^) 


70a>  0  5  if  0  <  oj  <  ujq 
UJ~2  if  U!  >  Uq  , 


(7.5) 


where  uq  =  0.007  (cycles  per  day).  For  continuity  of  Dn(k)  at  u> o,  70  =  u;0  1-5 . 
Equation  (7.5)  yields 


Dn-W(u))du)  —  429. 


(7.6) 


For  eastward-going  motions  (k  >  k0,u  >  0) 

Choosing  the  Bn~e(k)  to  fit  the  observed  zonal-wavenumber  spectrum  of  eastward¬ 
going  motions 


Bn-e(k) 


a2\k\  0  5  if  ko  <  k  <  k2 
|fc|~4  if  k  >  k2  , 


(7.7) 


where  k2  =  0.004  (cycles  per  km).  For  continuity  of  Bn(k)  at  k2)  a2  =  k2  3'5.  Equation 

(7.7)  gives 


=  2.7  x  107. 


(7.8) 
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Choosing  the  Dn-e( u)  to  fit  the  observed  frequency  spectrum  of  eastward- going 


motions 


Dn—e(s*j)  — 


7i  u  0  5  if  0  <  u)  <  cuq 
a;-1'5  if  u)  >  w0  , 


where  wq  =  0.007  (cycles  per  day).  For  continuity  of  Dn(k )  at  to0,  —  w0  1.  It  follows 
from  equation  (7.9)  that 


r+ °o 

/  Dn_e(u))du  =  48. 
Jo 


(7.10) 


For  standing  motions  (~k0  <  k  <  k0,co  >  0) 

Because  the  observations  provide  no  information  on  the  wavenumber  spectrum  at 
wavenumbers  —  k0  <  k  <  k0,  I  simply  assume  that  the  energy  of  standing  motions  is 
evenly  distributed  within  —k0<k<k0: 


Bn_0(k)  =  1  if  -k0  <  k  <  k0, 


(7.11) 


which  gives 


[ k°  Bn_0(k)dk  =  2k o  =  0.0008. 
J  —  kn 


(712) 


Choosing  the  Dn__o(uj)  to  fit  the  observed  frequency  spectrum  of  standing  motions 


Dn~  o(^) 


if  0  <  U  <U1\ 
if  u>  >  ui\  , 


(713) 


where  lo\  =  0.005  (cycles  per  day).  For  continuity  of  Dn(k)  at  a q,  72  =  u>1 125. 
Equation  (7.13)  gives 


f  +  OO 

/  Dn_0(o;)da;  =  1000. 

Jo 


(7.14) 


By  definition,  the  partitioning  of  the  variance  of  sea  surface  height  among  eastward- 
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going,  westward-going  and  standing  motions  for  the  model  is 


r+ oo  r- |-oo 

O rje  *  &r)W  ’  ~  Jk  Jq  Qn~e^n~e{^y)Dn—e(<(jJ^d(jjdk 

/—ko  /*+ oo  r+fco  T-foo 

I  Qn—w^n—w{'k'')Dn—w(^>)d'(jjdk  .  I  I  ^n—O-^Ti—O (^) -^n— 0 . (7. 1 5) 

-oo  */0  J—ko  JO 


Substitution  of  equation  (7.3)  to  (7.13)  into  equation  (7.15)  gives 


a%  :  a‘m  :  a';,,  =  1.5  x  109Q„_o  :  3.8  x  1010Q„_„  :  0.2Q„_C. 


(7.16) 


The  observations  in  chapter  4  show  that  the  energy  of  eastward-going,  westward-going 
and  standing  motions  is  approximately  in  the  ratio  of  1:2:4.  If  we  choose  Qn-e  —  1.0, 
then  according  to  equation  (7.16),  Qn-W  =  0.084  and  Qn-o  =  6.5  x  109. 

Choosing  the  Cn(l )  to  fit  the  observed  meridional-wavenumber  spectrum  of  sea 
surface  height 

fio  if  \l\  <  k 

Air 1/2  if/0<iii<ii 

r(o  =  r  (7-17) 

Air572 

i  r4  ^  i*i  >  *2 , 


(7-17) 


where  l0  =  0.0004;  /j  =  0.0025;  and  Z2  =  0.008  cycles  per  km.  For  continuity  of  Cn(l ) 
at  l\  and  /2 ,  A  =  *23//2*f2  a°d  A  =  /23^2-  Motions  within  |I|  <  /o  are  treated  as 
standing  waves  and  it  is  assumed  that  the  energy  is  constant  within  — 10  <  l  <  l0. 
Equation  (7.17)  gives 


f+OO 

/  Cn(l)dl  =  2.0  x  107, 

J  In 


(7.18) 


/  C„(l)dl  =  Pola  =  0.0004ft. 

JO 


(7.19) 


The  ratio  of  the  energy  of  propagating  motions,  which  includes  the  energy  of  northward- 
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going  and  southward-going  motions,  to  that  of  standing  motions  for  the  model  is 


afpT  :  alt  =  2.0  x  107  :  0.0004/?0.  (7.20) 

The  measurements  in  chapter  4  show  that  the  ratio  of  the  total  energy  of  northward¬ 
going  and  southward-going  motions  to  that  of  standing  motions  is  roughly  1:1.  Ac¬ 
cordingly, 


A)  =  5.0  x  1010. 


(7.21) 


Equation  (7.3)  to  (7.17)  yield 


/+ oo  r+oo  /•-foo 

/  /  $(fc,  Z,  a;,  n,  </>,  X)dkdldco  =  1.46  x  1O18J(0,  X)E0(n).  (7.22) 

-oo  J — oo  J —  OO 


The  energy  level  of  each  mode 

The  total  kinetic  energy  per  unit  surface  area  for  each  mode 

/-f-OO  Z'+OO  r+oo 

/  /  2tt2 (k2 +  l2)^(k,  l,  co,  n,(j),\)dkdldu; 

-oo  J ~ oo  J  —  oo 

=  3.0  x  I0ul((f>,  X)Eo(n).  (7.23) 


As  seen  in  chapter  4,  the  vertical  structure  of  eddy  kinetic  energy  can  be  simply 
represented  by  the  barotropic  and  first  two  baroclinic  modes  in  the  middle  of  the 
ocean  and  roughly  speaking,  the  kinetic  energy  of  the  first  three  modes  is  in  the  ratio 
of  l:l:l/2.  Accordingly,  if  we  choose  E’o(O)  =  1,  then  E'o(l)  =  1  and  E0(2)  =  0.5. 
Choose  the  I (cf),  A)  to  fit  the  observed  surface  kinetic  energy 
The  surface  kinetic  energy  for  the  model  is 

[  r+oo  /-foo  r+ oo  1 

Ek{(j),\,z  =  0)  =  T,Fn(z  =  0)\  /  /  2n2{k2 +  l2)<f>{k,l,oo,n,<P,\)dkdldco 

n—Q  L*/— OO  J  —  OO  J  —  OO 

=  9.9  x  10n/(^,A).  (7.24) 
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Equation  (7.24)  yields 


m  a) 


Ek(<l>,\z  =  0) 

9.9  x  1011 


(7.25) 


The  potential  energy  per  unit  surface  area  for  each  mode 


\  r+oo  r+oo  r+oo 

PE(n,  <f>,  A)  =  — —  /  /  /  $(kj,w,n,<f>,\)dhdldt 

ZltZ  J  —  00  J—  oo  J  —  oo 


1.5  x  1018 / ((j),X)E0(ji) 


(7.26) 


The  ratio  of  kinetic  energy  to  potential  energy  per  unit  surface  area  for  first  baro- 
clinic  mode 

If  we  take  /  =  10-4s-1,  R\  —  25  km  and 


A) 

PE(l,cf>,\) 


=  0.26. 


(7.27) 


The  variance  of  sea  surface  height 


^OM)  =  Y. 


PF*(z  =  0)  r+oo  r+ oo  r+oo 


/+oo  r+oo  r+oo 

/  /  $(k,l,ui,n,(t>,X)dkdldo 

-oo  J  —  OO  J—OO 


=  1.1  x  1012sm2^  /(</>,  A)  =  l.lsm20  2?fc(<^,  A,  z  =  0)  [cm2].  (7.28) 

The  zonal-wavenumber  spectrum  of  westward-going  motions  of  sea  surface  height 
(k  <  —k0,w  >  0) 


n—2  f2fp2(z  —  Q)  r+oo  r+oo 

Tw(A:,0,A)  =  - -  /  Qn-wEn-w{k)Dn^w(u>)Cn(l)Eo(n)I((j),  \)dldu 

“n  Q  JO  J- OO 


=  2.1  x  103sm2^>  /(</>,  \)Bn-w(k)  [cm2/(cycles/km)j. 


(7.29) 


The  zonal-wavenumber  spectrum  of  eastward- going  motions  of  sea  surface  height 
(k  >  k0,w  >  0) 


11+  f2F2(z  =  0)  /’+00  r+oo 
T„e((c>,A)  =  £i_»L - 1  f  1 

ZA  9  •'»  7-00 


n+oo 

Qn-eBn-e(k)Dn-e(u)Cn(l)E0{n)I(<l>,\)dldu 

-OO 


=  2.8  x  103 sin2 <f)  I ((f,\)Bn-e(k)  [cm2/(cycles/km)]. 


(7.30) 


The  scalar  zonal- wavenumber  spectrum  of  sea  surface  height  ( k  >  k0) 


/+OO  /"foo 

/  FT)(k,l,u,(f),\)dldu} 

-OO  J  —  OO 

B-qu; ( k ,  A)  +  (f>,  A) 

=  sin2  <f>  7(0,  A)  [2100i?n_ti,(A;)  +  2800£?7l_e(fc)]  [cm2/(cycles/km)].  (7.31) 

The  scalar  meridional-wavenumber  spectrum  of  sea  surface  height  ( l  >  l0) 

TEz2  f2F2(z  =  O')  r+°o  r+ oo 

T^/.^A)  =  - -  /  nn(fc,£4;)Cn(/)£?o(w)/(0,  X)dldco 

n=0  9  J~°°  J- oo 

=  1.3  x  104 sin2 0  7(0,  X)Cn(l)  [cm2/(cycles/km)].  (7.32) 


The  frequency  spectrum  of  westward-going  motions  of  sea  surface  height  (k  < 
—k0,uj  >  0) 


T A) 


f2F2(z  =  0)  f+°°  r~k o 

- L  /  Qn_wBn^w(k)Dn_w(u))Cn(l)Eo(n)I (f>,  X)dkdl 

71=0  9  J~°°  ^-OO 

3.5  x  108  sin2  (f>  \)Dn-w(u>)  [cm2/(cycles/day)].  (7.33) 


The  frequency  spectrum  of  eastward- going  motions  of  sea  surface  height  (k  > 
k0,u  >  0) 

n=2  f2jT2(z  —  0)  r+oo  r+o o 

T^e(o;,  0,  A)  =  - -  /  Qn-eB„-e(k)Dn-e(w)Cn(l)E0(n)I(f),  \)dkdl 

n=o  9  J-°°  Jko 

=  1.6  x  109 sin2 <f>  I(f>,  \)Dn-e(u)  [cm2/(cycles/day)j.  (7-34) 


The  frequency  spectrum  of  standing  motions  of  sea  surface  height  (—k0  <  k  < 
k0,u>  >  0) 


71=0  9 

3.0  x  10 ssin2(j)  I(<t>,\)Dn_0(w) 


/+oo  rko 

/  Qn-oBn-o{k)Dn_0(u)Cn(l)E0(n)I((f),  X)dkdl 

-OO  V—  &0 

[cm2/(cycles/day)l.  (7.35) 


The  scalar  frequency  spectrum  of  sea  surface  height  (w  >  0) 


^  T)W  (oj,  (j>,  A)  +  T jje(u),  </>,  A)  +  T^o (c+  A) 

=  sin2cf>  /(</>,  A)  [3.5  x  108Dn_w (to)  +  1.6  x  109£>„_e(u;)  +  3.0  x  1087An_0(w)]  ,(7.36) 

in  the  units  of  cm2/(cycles/day). 

The  scalar  zonal-wavenumber  spectrum  of  temperature  ( k  >  ko ) 

n=2  r)f)n  r+oo  r+oo 

T*(M,<M)  =  Y,f2Gl(z)(^-)2Eo(n)I((t>,X)  /  /  U(k,u)Cn(l)dldto 

n_Q  uZ  J— oo  J —oo 

=  6.1  x  10 3sin2(f)  7(0,  A)(^)2  [Gl(z)  +  0.5  x  G22{zj\ 

[Bn~w(k)  +  1.4  x  Bn.e(k )]  [°C2/ (cycles/km)].  (7.37) 

T/*e  scalar  meridional-wavenumber  spectrum  of  temperature  (l  >  lo ) 

^  Z’+OO  /*+oo 

Fe(l,z,  <f>,  A)  =  E/2^)^)2  /  /  n(fc,u;)dfcdu;C'n(0^oH/(^A) 

n=0  ^  _°°  ^-oo 

no 

=  3.8  x  104sm2</>  7(^>,  A)(-7p-)2[G2(z)  +0.5  x  G2(2;)]C'„(/)  [°C2/(cycles/km)].  (7.38) 

1//6 

77ie  frequency  spectrum  of  zonal  velocity  (to  >  0) 

/’H-OO  r+OC 

Ttt(u/,  z,<f>,  A)  =  E  /  /  47r2i2F2(z)$(7,  /,  w,  n,  <j>,  \)dkdl 

n= 0J-°°  J~°° 

=  1.25  x  10n/(<£,  A)  [f2(^)  +  F2(z)  +  0.5 Ff(z)] 

x  [Z>„_w (to)  +  4.5  x  Dn^e(to)  +  0.87  x  Dn^(to)}  [cm2/s2/ (cycles/day)](7.39) 
The  frequency  spectrum  of  meridional  velocity  (to  >  0) 

rc— 2  r+oo  /•+ oo 

T„(a;,  z,  (f),  A)  =  E  /  /  4TT2k2F2(z)$(k,l,to,n,(j),  \)dkdl 

n-  0  7-oo 

=  1.3  x  10n/(</>,  A)  [F02(^)  +  F2(z)  +  0.5-F2(z)] 
x  [D„_*(a;)  +  8.7  x  £>n_e(a;)  +  7.0  x  IO^TVoM]  ,  (7.40) 
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in  the  units  of  cm2/s2/(cycles/day). 

The  frequency  spectrum  of  temperature  (co  >  0) 

A)  =  ^  f2Gn(z)(-~ -)2  f  f  $(k,l,u),n,(f>,\)dkdl 

n=Q  OZ  J  ~  oo  J  ~  oo 

=  1-0  X  109 sin2 f)  I{4 >,  A)[G2(z)  +  0.5  x  G2{z)}{^)2 
x  [ Dn-W(u >)  +  4.5  x  Dn-e{ u)  -f  0.85  x  Dn^0(u>)]  [°C2/(cycles/day)].  (7.41) 

7.2  Model  and  data  comparison 

As  in  chapter  5,  the  model  spectra  will  be  compared  with  the  corresponding  obser¬ 
vations.  The  data  used  here  are  the  same  as  those  used  in  chapter  5. 

The  frequency /wavenumber  spectrum  of  sea  surface  height 

The  frequency  and  wavenumber  spectra  of  westward-going  and  eastward-going 
motions  in  sea  surface  height  are  displayed  in  figure  7-1.  The  model  spectra  fit  the 
observations  quite  well  in  terms  of  both  spectral  shape  and  energy  level  (figure  7- 
1).  As  shown  in  figures  7-2  and  7-3,  the  scalar  zonal- wavenumber  spectrum  and 
meridional-wavenumber  spectrum  of  the  model  are  in  rough  agreement  with  the  cor¬ 
responding  observations.  The  energy  level  of  model  frequency  spectra  in  figures  7-2 
to  7-5  are  systematically  higher  than  the  observations  and  the  differences  are  within 
a  factor  of  2.  Figure  7-4  shows  that  superimposed  on  the  background  continuum, 
there  is  a  striking  peak  at  the  annual  period  on  the  observed  spectrum.  The  peak 
associated  with  annual  cycle  only  appears  on  the  frequency  spectrum  of  sea  surface 
height:  there  is  no  significant  peak  at  the  period  of  one  year  on  the  frequency  spec¬ 
tra  of  horizontal  velocities  and  temperature.  The  seasonal  cycles  should  be  treated 
exclusively.  This  is  left  for  future  study. 

The  scalar  wavenumber  spectrum  of  temperature 

The  scalar  zonal-wavenumber  spectrum  of  temperature  is  shown  in  figure  7-6.  As 
shown  there,  the  agreement  between  the  model  spectrum  and  the  observed  spectrum 
is  quite  pleasing. 
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The  scalar  frequency  spectrum  of  horizontal  velocities  and  temperature 
Figure  7-7  displays  the  frequency  spectra  of  two  components  of  horizontal  veloc¬ 
ities  and  temperature  at  5  different  depths.  Because  the  data  sets  are  short,  the 
spectra  have  large  confidence  intervals.  The  model  spectra  are  indistinguishable  from 
the  observations  with  respect  to  a  95%  confidence  interval. 


7.3  Summary 

The  directional  form  for  the  energy  density  $(k,  l,  u,  n ,  0,  A)  is  produced  in  this  chap¬ 
ter.  Compared  with  the  model  produced  in  chapter  5,  this  one  can  differentiate 
westward-going  energy  from  eastward-going  energy;  however  it  is  more  complex.  Only 
the  first  three  modes  are  included  in  my  model:  the  barotropic  and  first  two  baroclinic 
modes.  For  each  mode  I  choose 


$(&,  l,  u),  n,  0,  A)  =  fl n(k,  u)Cn{l)E0(n)I(f>,  A). 


(7.42) 


fIn(A;,  oj)  has  distinct  forms  for  eastward-going,  westward-going  and  standing  motions, 


Tin  (&?  kj)  —  ^ 


Qn— wEn— w{T)Dn— w  (ca) 


Q 


n—eBn—eih')T)n—e(uf) 


Qn—oEn—o  (k)Dn_0  (ca) 


for  k  <  —k0,  u  >  0 
for  k  >  k0,u>  >  0 
for  —  ko  <  k  <  ko,  u>  >  0 


(7.43) 


where  k0  =  0.0004  cycles/km.  Based  on  various  observations,  I  obtain 
(1)  the  zonal- wavenumber/frequency  spectral  shape  IIn(fc,a;) 

For  westward-going  motions  ( k  <  —k0,u  >  0) 


Bn—w  (^)  —  \ 


8.6  x  108|A:|-0-5 

\k\~4 


if  —ki  <  k  <  —ko 
if  k  <  k\  , 


(7.44) 
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where  k\  =  0.0028  cycles  per  km. 


Dji— tu(^)  — 


1.7xl0^-°'5  if  0  <  w  <  o;0 
a;-2  if  a:  >  u0  , 


(7.45) 


where  u>q  —  0.007  cycles  per  day. 


Qn-w  =  0.084 


(7.46) 


For  eastward- going  motions  ( k  >  k0,oo  >  0) 


Bn-e(k)  = 


2.5  x  108|&|~°'5  ifk0<k<k2 
|ft|-4  if  k  >  k2  , 


(7.47) 


where  k2  —  0.004  cycles  per  km. 


143a;  0  5  if  0  <  ui  <  cv0 
a j~15  if  uj  >  too  , 


(7.48) 


where  ujq  =  0.007  cycles  per  day. 


Qn-e  =  1-0 


(7.49) 


For  standing  motions  (— A:0  <  k  <  k0,  uj  >  0) 


Bn_0(k )  =  1  if  —k0  <  k  <  k0 


(7.50) 


Dn-o(w)  — 


where  u>i  =  0.005  cycles  per  day. 


752a;-0-75  if  0  <  u  < 
uj~2  if  uj  >  uii  , 


(7.51) 


Qn-0  =  6.5  X  109. 


(7.52) 
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(2)  the  meridional-wavenumber  spectral  shape 


’  2.5  x  1011  if  |/j  <  lQ 

I  2.2  x  108!/!-1/2  if  l0  <  \l\  <  h 

1.4  x  103|/|— 5/2  if  i,  <  |/|  <  l2 

.  r4  « w  >  k , 

where  l0  =  O.OOOl,  l\  =  0.0025,  and  I2  =  0.008  cycles  per  km. 
(3)  the  partition  among  vertical  modes 

£0(0)  -  1,  J50(l)  =  landE0{2)  =  0.5. 


(7.53) 


(7.54) 


where  A  is  the  longitude  ,  (f>  is  the  latitude  and  the  units  of  Ek  are  cm2/s2. 

The  scalar  formula  for  &(k,l,u,n,  <j>,  A),  which  is  produced  in  chapter  5,  can 
be  obtained  from  the  directional  form  by  adding  up  the  energy  of  eastward-going, 
westward-going,  and  standing  motions. 
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Figure  7-1:  The  zonal-wavenumber  spectrum  of  westward-going  motions  (a)  and 
eastward-going  motions  (b),  and  the  frequency  spectrum  of  westward- going  motions 
(c)  and  eastward-going  motions  (d)  of  sea  surface  height.  The  spectra  of  measure¬ 
ments  and  model  are  plotted  in  solid  lines  and  dotted  lines,  respectively.  The  peak 
at  the  period  about  60  days  in  figure  (c)  and  (d)is  due  to  tidal  aliasing. 
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west-east  wavenumber  (cydes/km) 


Figure  7-2:  Scalar  zonal-wavenumber  spectrum  of  sea  surface  height.  Solid  line  is  the 
observed  spectrum.  Dotted  line  is  the  corresponding  model  spectrum. 


north-south  wavenumber  (cyctes/km) 

Figure  7-3:  Scalar  meridional-wavenumber  spectrum  of  sea  surface  height.  Solid  line 
is  the  observed  spectrum.  Dotted  line  is  the  corresponding  model  spectrum. 


Figure  7-4:  Frequency  spectrum  of  standing  motions  of  sea  surface  height.  Solid 
line  is  the  observed  spectrum.  Dotted  line  is  the  corresponding  model  spectrum. 
Superimposed  on  the  background  continuum,  there  is  a  striking  peak  at  the  annual 
period. 


Figure  7-5:  Scalar  frequency  spectrum  of  sea  surface  height.  Solid  line  is  the  observed 
spectrum.  Dotted  line  is  the  corresponding  model  spectrum.  The  observed  spectrum 
is  dominated  by  the  annual  cycle. 
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Figure  7-7:  Scalar  frequency  spectra  of  zonal  (a),  meridional  velocity  (b)  and  tem 
perature  (c)  at  depths  of  160m  (i),  580m  (ii),  595m  (iii),  1230m  (iv)  and  3000m  (v) 
The  solid  line  is  the  observed  spectrum  and  the  dotted  line  is  the  model  spectrum. 


Chapter  8 


The  observed  dispersion 
relationship  for  North  Pacific 
Rossby  wave  motions 

8.1  Introduction 

Most  of  this  chapter  was  published  by  Zang  and  Wunsch  (1999).  In  the  preceding 
chapters,  efforts  are  made  to  produce  a  universal  continuous  spectrum  for  low  fre¬ 
quency  variability.  The  low  frequency  variability  consists  of  Rossby  waves,  eddies  and 
etc.  Rossby  waves  are  important  parts  among  the  continuous  spectrum  we  produced 
in  previous  chapters.  What  percentage  of  the  ocean  variability  can  be  explained  by 
the  linear  Rossby  wave  theory  is  of  interest  not  only  for  their  links  to  the  forecast  of 
low  frequency  oceanic  variability  but  also  from  a  theoretical  point  of  view.  Rossby 
waves  play  an  important  role  in  maintaining  the  western  boundary  layer  and  the 
transient  adjustment  of  ocean  circulation  to  changes  in  atmospheric  forcing.  The 
presence  of  the  first  baroclinic  long  Rossby  waves  was  found  early  from  upper-ocean 
temperature  data  (White,  1977;  Kessler,  1990).  White  (1977)  analyzed  the  XBT  data 
for  the  tropical  North  Pacific  and  discovered  that  the  annual  vertical  displacements 
in  the  main  thermocline  propagate  westward  as  long  baroclinic  waves,  with  a  phase 
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speed  twice  that  of  free  first  baroclinic  Rossby  waves.  Recently,  Chelton  and  Schlax 
(1996)  performed  an  analysis  of  the  3-year  high  quality  TOPEX/POSEIDON  altime¬ 
ter  data.  They  filtered  out  high  frequency  and  small  scale  variabilities  to  highlight 
the  long  Rossby  waves  and  found  that  the  disagreement  between  the  apparent  phase 
speeds  of  Rossby  waves  and  the  theoretically  predicted  values  is  universal.  Poleward 
of  the  tropical  band  from  10°S  to  10°N,  the  apparent  phase  speeds  are  systematically 
greater  than  those  predicted  by  standard  theory  for  free  first  baroclinic  long  Rossby 
waves  (Figure  8-1). 

The  difference  between  the  apparent  phase  speeds  and  those  predicted  by  standard 
theory  has  been  attributed  by  several  authors  (e.g.,  Qiu  et  ah,  1997;  Killworth,  1997) 
to  the  breakdown  of  one  or  more  of  the  assumptions  inherent  in  the  linear  theory. 
Possible  mechanisms  include  atmospheric  forcing,  wave  and  mean  flow  interactions, 
wave  and  wave  interaction,  coupling  of  vertical  modes  and  topography.  Some  of  these 
factors  will  be  examined  in  the  following. 

White  (1977)  studied  the  ocean  response  to  annual  wind  forcing.  He  showed  that 
when  superimposed  on  the  local  response,  the  free  long  Rossby  wave  yields  a  long 
wave  which  propagates  at  a  phase  speed  twice  that  of  the  free  long  Rossby  wave. 
Qiu  et  al.  (1997)  extended  the  theory  by  White  (1977).  They  investigated  what 
the  phase  speed  of  the  total  response  was  when  in  the  presence  of  the  forced  wave 
response  to  wind  forcing  in  the  interior  ocean,  the  eastern  boundary  generated  free 
Rossby  wave  and  dissipation.  They  showed  that  the  forced  wave  pattern  propagates 
westward  at  2cr,  where  cr  is  the  speed  of  the  free  first  baroclinic  long  Rossby  wave.  In 
low  latitudes  the  westward  phase  speed  of  the  total  response  favors  the  speed  of  lcr 
because  of  the  boundary-generated  free  Rossby  waves.  Toward  higher  latitudes,  the 
boundary-generated  free  waves  become  highly  dissipative,  the  interior-forced  wave 
motion  becomes  dominant.  Thus  the  total  response  tends  to  propagate  at  2cr. 

Killworth  et  al.  (1997)  investigated  the  effect  of  mean  zonal  flow  on  the  propaga¬ 
tion  of  first  baroclinic  long  Rossby  waves.  Using  historical  hydrographic  data,  they 
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obtained  the  climatological  buoyancy  frequency  profiles  and  the  vertical  shear  of  the 
mean  flow.  They  found  that  the  effect  of  the  mean  barotropic  zonal  flow  is  purely  a 
Doppler  shift  and  the  inclusion  of  the  barotropic  flow  has  small  and  mixed  effects  on 
the  phase  speed  of  the  first  baroclinic  Rossby  wave.  The  major  changes  to  the  Rossby 
wave  phase  speed  are  caused  by  the  presence  of  second  baroclinic  zonal  mean  flow, 
which  modifies  the  potential  vorticity  gradient.  Assuming  some  idealized  form  for 
the  model  buoyancy  and  mean  shear  flow,  Dewar  (1998)  performed  a  similar  study 
and  found  that  baroclinic  waves  are  accelerated  through  the  interactions  of  the  waves 
with  the  background  potential  vorticity  field,  de  Szoeke  and  Chelton  (1999)  further 
investigated  what  features  in  the  mean  flow  and  stratification  cause  the  acceleration 
of  the  long  baroclinic  Rossby  waves.  By  studying  a  simple  three-layer  model,  they 
concluded  that  the  modification  of  long  baroclinic  Rossby  waves  is  associated  with 
potential  vorticity  homogenization. 

So  far,  all  studies  attribute  the  difference  between  observations  and  theory  to  the 
limitation  of  the  simple  linear  theory.  One  important  fact  has  been  neglected.  The 
observed  phase  speed  obtained  by  Chelton  and  Schlax  (1996)  is  not  the  phase  speed 
of  free  first  baroclinic  Rossby  waves.  As  seen  in  chapter  4,  the  sea  surface  height 
variability  is  a  broad  band  process.  It  is  made  up  of  many  different  signals  such 
as  free  Rossby  waves,  forced  waves,  etc.  Moreover,  although  in  general  sea  surface 
height  variability  mainly  reflects  the  first  baroclinic  mode,  the  barotropic  mode  and 
high  baroclinic  modes  play  an  important  role  over  many  parts  of  the  ocean  (Wun- 
sch,  1997).  Therefore,  not  all  the  sea  surface  height  variability  is  necessarily  first 
baroclinic  Rossby  waves.  The  observed  phase  speed  obtained  by  Chelton  and  Schlax 
(1996)  is  not  the  phase  speed  of  free  first  baroclinic  Rossby  waves,  but  the  phase 
speed  of  the  total  response  which  is  composed  of  barotropic  Rossby  waves,  first  and 
second  baroclinic  Rossby  waves,  forced  waves,  etc.  Although  the  linear  Rossby  wave 
theory  is  a  first-order  good  approximation  to  oceanic  low  frequency  variability,  the 
actual  ocean  is  turbulent  and  much  of  the  actual  time-dependent  motion  can  not 
be  adequately  described  by  linear  wave  models.  How  to  separate  the  free  Rossby 
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waves  from  other  waves  and  to  which  extent  the  linear  Rossby  wave  theory  is  a  good 
representation  of  the  true  ocean  variability  are  the  incentive  for  this  work.  Because 
the  standard  theory  for  free  Rossby  waves  is  expressed  in  terms  of  frequency  and 
wavenumber,  the  best  way  to  compare  observations  with  theory  is  to  express  the 
data  in  Fourier  space,  that  is,  analyze  the  three  dimensional  spectrum  of  the  data. 
The  determination  of  the  frequency/wavenumber  structure  of  the  wave  field  is  one 
of  the  most  powerful  techniques  for  the  study  of  the  properties  of  ocean  variability. 
The  frequency/wavenumber  spectrum  yields  the  direction  and  phase  speed  of  wave 
propagation,  and  therefore  helps  us  to  forecast  the  propagating  signals  and  clarify 
the  fundamental  processes  of  wave  generation.  It  can  also  help  us  separate  different 
physical  processes  by  taking  advantage  of  the  fact  that  different  waves  have  different 
frequency  and  wavenumber  characteristics.  Cipollini  et  al.  (1997)  performed  Fourier 
transforms  of  TOPEX/POSEIDON  altimeter  data  and  ERS-1  along-track  scanning 
radiometer  sea  surface  temperature  data  near  34°N  in  the  Northeast  Atlantic.  They 
observed  westward  propagating  signals  with  three  distinct  speeds.  However,  they 
neither  compared  the  observations  with  the  theory  nor  provided  estimates  of  uncern- 
tainty.  Using  a  filter  in  time  and  a  wavelet  transform  in  space,  Wang  et  al.  (1998) 
analyzed  the  gridded  sea  surface  height  data  and  identified  two  types  of  Rossby  waves 
with  different  frequencies  and  wavenumbers  along  25° S  in  the  South  Pacific.  No  esti¬ 
mates  of  uncertainty  was  provided.  Here  we  apply  array  signal  processing  technique 
to  the  original  along  track  sea  surface  height  data.  We  will  provide  estimates  of  un¬ 
certainty  and  compare  the  observations  with  the  theoretical  dispersion  relationships 
for  Rossby  waves. 

The  presentation  of  this  chapter  is  organized  as  follows.  The  array  processing 
method  will  be  reviewed  in  section  2.  Section  3  presents  data  and  data  processing. 
The  results  are  discussed  in  section  4.  A  brief  summary  is  drawn  in  section  5. 
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Figure  8-1:  (A)  Gloabally  distributed  estimates  of  the  phase  speeds  of  westward- 
propagating  sea  level  signals.  The  solid  circles  correspond  to  Pacific  estimates,  and 
the  open  circles  correspond  to  Atlantic  and  Indian  Ocean  estimates.  The  continuous 
line  stands  for  the  global  average  latitudinal  variation  of  the  phase  speed  predicted  by 
the  standard  theory  for  extratropical  freely  propagating,  linear,  first  baroclinic  mode 
long  Rossby  waves.  (B)  Ratio  of  the  observed  phase  speeds  to  the  phase  speeds  pre¬ 
dicted  by  the  standard  theory  at  the  same  geographical  locations  as  the  observations 
(Chelton  and  Schlax,  1996). 
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8.2  Beam- forming  method 

Two  types  of  beam-forming  methods  will  be  used  here,  that  is,  classic  beam-forming 
technique  and  minimum  variance  estimator.  Classic  beam-forming  techniques  are  a 
“delay- and-sum”  form  of  processing  in  the  Fourier  frequency/wavenumber  domain. 
Data,  Xi  ( t ) ,  at  each  position  are  temporally  band-pass  filtered,  by  ordinary  fast- 
Fourier  transform,  in  a  narrow-band  surrounding  temporal  frequency  a ,  generating  the 
Fourier  transforms  Xi  (to) .  All  cross-power  combinations  Rij  (tu)  =<  Xi  (u)xj  (u>)*  > 
are  formed  (the  *  denotes  complex  conjugation  and  the  bracket  average  is  estimated 
by  frequency-band  averaging).  Denote  the  vector  separation  between  positions  i,j  as 
Arjj.  The  estimated  dominant  wavenumber  present  is  obtained  from  the  maximum 
over  all  plausible  wavenumbers,  k  of 

$  (uj,  k)  =  $  (u,  k,l)  =  ^2  WiWjRij  (to)  exp  (-ik  •  Ar^)  ,  (8.1) 

ij 

where  Wi  are  weights  which  are  used  to  control  the  shape  of  the  wavenumber  win¬ 
dow.  The  “beam-pattern”  is  the  field  $  over  k  evaluated  in  the  presence  of  a  single, 
monochromatic  plane  wave  exp  (ik0  •  r  —  iuit) .  Unlike  methods  based  upon  multi¬ 
dimensional  fast  Fourier  methods,  the  spatial  separations  ArZJ  can  be  quite  irregular. 
As  demonstrated  by  Capon  (1969),  equation  (8.1)  is  a  biased  estimator  of  the  actual 
three-dimensional  spectrum  as  a  result  of  window  tapering.  Therefore,  corrections 
must  be  made  to  compensate  for  the  loss  of  energy  due  to  window  tapering  (see 
Capon,  1969).  Oceanographic  applications  of  this  and  related  methods  may  be  seen, 
e.g.  in  Wunsch  and  Hendry  (1972),  and  elsewhere. 

The  basic  frequency  resolution  of  this  form  of  analysis  is  given  by  the  reciprocal  of 
the  data  duration,  reduced  by  the  frequency-band  averaging  as  in  an  ordinary  power 
density  spectral  estimate;  the  highest  frequency  estimated  is  given  by  the  normal 
Nyquist  criterion  (data  are  uniformly  spaced  in  time).  Wavenumber  resolution  is 
similarly  determined  by  the  zonal  array  dimension  with  the  spatial  Nyquist  deter- 
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mined,  theoretically,  by  the  smallest  available  zonal  spatial  separation.  Because  the 
data  positions  need  not  be  uniformly  spaced,  the  actual  practical  Nyquist  wavenum¬ 
ber  is  dependent  upon  the  array  geometry  and  signal-to-noise  ratio. 

The  minimum  variance  estimator  is  a  high  resolution  array  processing  method 
which  is  due  to  Capon  (1969).  It  is  sometimes  termed  the  “maximum-likelihood 
method”,  because  the  form  of  this  estimator  is  similar  to  that  in  the  maximum- 
likelihood  estimation  of  the  amplitude  of  a  sine  wave  of  known  frequency  in  Gaussian 
random  noise  (Capon  et  al.,  1967).  The  array  output  power  contains  not  only  the 
desired  energy  along  the  look  direction  but  also  the  undesired  one  from  other  direc¬ 
tions.  The  minimum  variance  estimator  adaptively  constructs  a  weighting  function  to 
minimize  the  array  output  power  while  maintaining  the  gain  along  the  look  direction 
to  be  unity  (Pillai,  1989). 

Let  R  be  the  cross-power  spectral  matrix  and  w  be  the  weighting  vector  to  be 
specified,  then  the  output  power  is 

S(u)  =  w+Rw.  (8.2) 

Capon’s  Estimator  is  equivalent  to  the  following  problem  (Pillai,  1989) 

min  w+Rw  subject  to  w+a  =  1  (8-3) 

where  a  is  the  steering  vector 


a=  [eifel,eifel,...eifeN]T 


(8.4) 


The  constrained  minimization  can  be  found  by  using  a  Lagrangian  multiplier  A  to 
construct  an  new  function 


T  =  w+Rw  +  A(w+a  —  1)  '  (8-5) 
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If  the  matrix  R  is  positive-definite,  the  optimum  weight  vector  is  readily  given  by 


R  *a 
a+R-V 


(8,6) 


The  corresponding  output  power  with  this  optimum  weight  vector  is 


$(u;,  k,  l)  = 


1 


1 


a+Qa 


N 

E  Qji{u)  exp [ik(xj  -  £;)] 

j,i= i 


(8.7) 


which  is  defined  as  the  three  dimensional  estimator.  Q  is  the  inverse  of  the  matrix 
R.  Assuming  the  sensor  outputs  to  be  stationary  Gaussian  processes,  Capon  and 
Goodman(1970)  showed  that  the  random  variable  $(u ),k,l)  is  a  multiple  of  a  chi- 
square  variable  with  2 M  degrees-of-freedom,  where  M  is  the  number  of  frequency 
bands  averaged.  In  this  study,  we  choose  m  =  3  then  there  are  thus  about  six 
degrees-of-freedom. 

The  reader  beware  that  the  necessary  condition  for  equation  (8.7)  to  hold  is  that 
the  cross  spectral  matrix  R  must  be  positive-definite.  This  necessary  condition  is 
often  unattainable.  The  singular  covariance  matrix  is  a  well  known  problem  in  the 
array  signal  processing  literature  and  can  be  fixed  in  a  number  of  ways.  Diagonal 
loading  is  one  of  the  common  methods  (Baggeroer  and  Cox,  1999).  However,  when 
diagonal  loading  is  applied,  bias  will  be  introduced  into  the  spectral  estimates.  This 
issue  has  been  examined  recently  by  Baggeroer  and  Cox  (1999).  When  diagonal 
loading  is  applied  the  cross  spectral  matrix  R  is  replaced  by 


R(<5)  =  R  +  51, 


(8.8) 


where  5  represents  the  noise  power  which  is  added  to  make  the  spectral  matrix  non¬ 
singular.  With  this  cross  spectral  matrix,  the  weight  vector  of  the  minimum  variance 
beamformer  becomes 


w((5) 


R(<5)_1a 

a+R((5)-1a 


(8.9) 
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If  the  above  weight  vector  is  applied  to  the  unloaded  sample  cross  spectral  matrix, 
the  beam  output  is 

$i(u,k,l)  =  w(5)+Rw((5).  (8.10) 

On  the  other  hand,  if  the  loaded  cross  spectral  matrix  is  used,  the  beam  output  is 

$2(w,M)=  [a+R^-'a]-1.  (8.11) 

These  formulas  were  derived  by  Baggeroer  and  Cox  (1999).  There  are  no  analytic 
solutions  of  the  bias  due  to  diagonal  loading.  However,  the  bias  can  be  estimated 
through  Monte  Carlo  methods  (Baggeroer  and  Cox,  1999).  The  effect  of  diagonal 
loading  on  the  confidence  limit  of  minimum  variance  estimator  is  still  unclear  and  is  an 
outstanding  problem  in  array  signal  processing  (Baggeroer,  personal  communication). 

8.3  Data  and  data  processing 

Here,  following  CS,  we  use  TOPEX/POSEIDON(T/P)  altimeter  data,  from  the  pe¬ 
riod  1  January  1993  to  21  September  1997,  edited  and  corrected  as  described  by  King 
et  al.  (1994)  except  for  the  use  of  a  more  recent  tidal  model.  The  data  set  available 
to  CS  spanned  about  three  years.  Here  the  five  years  of  data  improves  both  the 
frequency  resolution  and  the  statistical  reliability.  T/P  provides  global  sea  surface 
height  observations  every  9.91  days.  Data  within  an  orbit  cycle  were  assigned  to  a 
common  time.  This  assignment  will  have  a  distorting  effect  on  the  result,  but  as 
will  be  seen,  the  results  suggest  it  is  not  of  first-order  importance.  The  time-mean 
was  removed  from  each  along-track  point  to  obtain  the  time  dependent  part,  then 
a  three  block  non-overlapping  time  average  was  performed  to  filter  out  all  motions 
with  periods  shorter  than  59.5  days.  The  result  is  an  array  of  elevation  anomalies 
rf  (0i,  Aj,t)  where  <j>  is  the  latitude,  A  the  longitude  and  t  =  lAt,  2Ai, ...,  58At,  is  a 
time  index,  with  At  =29.7  days.  A  large,  and  large-scale,  steric  component  exists  at 
the  annual  period  in  the  altimetric  data  (see  Stammer,  1997b);  because  dynamical 
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effects  at  one  year  are  of  interest  here,  this  contribution  was  not  removed,  and  we 
rely  upon  the  wavenumber  rejection  capability  of  our  method  to  separate  the  k  =  0 
contributions  from  others. 

One  of  the  difficulties  with  the  study  of  Rossby  waves  is  that  the  propagation  is 
not  only  anisotropic,  but  also  inhomogeneous.  In  particular,  because  the  meridional 
variation  of  f  (y)  is  not  linear,  and  i?„,  is  a  function  of  y  (and  x)  the  phase  velocities 
are  functions  of  latitude.  Phase  lines  initially  aligned  with  a  meridian  will  “turn” 
because  the  phase  velocity  is  greater  at  low  latitudes  than  at  high  latitudes  (see  Schopf 
et  al.,  1981).  If  the  ocean  occupied  the  entire  globe,  the  proper  basis  set  would  be 
the  spherical  harmonics.  In  the  present  situation  we  represent  the  waves  as  sines  and 
cosines  in  narrow  latitude  bands  in  a  series  of  boxes  depicted  in  Figure  8-2  and  listed 
in  Table  8-1.  The  latitude  range  is  restricted  to  0.2°  so  that  the  meridional  changes 
are  negligible.  The  longitude  bands  are  sufficiently  restricted  that  zonal  variations  in 
Rn  are  slight. 

Seven  boxes,  chosen  somewhat  arbitrarily  to  represent  low,  middle  and  high  Pacific 
latitudes  as  well  as  some  degree  of  east-west  dependence,  are  used.  Consider  any  box 
centered  at  latitude  0O-  Then  within  the  box,  we  can  think  of  each  position  of  an 
altimetric  time  series  as  though  it  were  part  of  an  array  (see  Fig.  8-2), 

{xi  (t)}  =  {rf  (fa,  X i,  t)}  ,  1  <  %  <  N,  1  At  <  t  <  58A t,  (8.12) 

where  Xi  ( t )  is  the  time  series  at  pseudo-sensor,  i,  and  analyze  it  by  standard  “beam¬ 
forming”  techniques  (e.g.  Capon,  1969;  Lacoss,  1971;  Bath,  1974).  An  alternative 
procedure  would  be  direct  three  dimensional  least-squares  fitting  of  sines  and  cosines; 
array  processing  procedures  are  however,  efficient  and  convenient. 

Because  of  the  “diamond-pattern”  in  which  the  satellite  ground-tracks  cover  the 
ocean,  we  can  choose  analysis  positions  <f>i,  Xi  within  the  boxes  either  with  a  regular 
or  irregular  spacing  along  and  between  these  tracks.  The  zonal  spacing  between 
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sensors  is  vanishingly  small  near  where  the  descending  and  ascending  arcs  intersect. 
So-called  data-adaptive  techniques  for  array  processing  are  well  known,  especially 
Capon’s  (1969)  method.  These  methods  are  most  powerful  when  arrays  are  small 
and  irregular,  and  the  plane-wave  signal-to-noise  ratio  is  large.  In  the  present  case, 
we  arrange  for  the  data  coverage  to  be  fairly  uniform  and  large-scale,  and  the  Capon 
method  did  not  produce  better  results  than  did  the  more  robust  and  conventional 
beamforming  method.  An  attempt  to  examine  high  wavenumbers  by  employing  a 
very  disparate  array  spacing  and  the  Capon  method  will  be  described  briefly  at  the 
end. 

After  filtering,  each  “sensor”  has  an  altimetric  time  series,  Xi(t),  at  intervals, 
At  =  29.7  days,  with  duration  T  =  58A t.  To  stabilize  the  cross-power  estimates, 
frequency-band  averaging  was  done  over  three  adjacent  values  which  thus  sets  the 
approximate  frequency  resolution  of  the  analysis.  Let  the  zonal  width  of  the  array 
(its  aperture)  be  L,  and  let  the  smallest  zonal  spacing  between  any  two  elements  be 
Ar  (in  kilometers).  Beamforming  produces  estimates  of  the  power  density,  $  (a,  k,  l ) , 
in  a  band  of  zonal  wavenumbers,  2n/L  <  k  <  7r/(A r).  Each  box  extends  sufficiently 
far  in  the  zonal  direction  to  distinguish  the  first  mode  baroclinic  wave  from  motions 
at  or  near  k  =  0.  This  wavenumber  band  is  separately  computed  for  each  of  the 
resolved  frequencies.  Because  the  meridional  aperture  is  so  small,  there  is  almost 
no  resolving  power  in  the  meridional  direction  (wavenumber  l)  which  is  therefore, 
indistinguishable  from  l  =  0  in  the  results. 
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Figure  8-2:  Rectangles  are  the  seven  regions  chosen  to  study  the  frequency  wavenum¬ 
ber  spectrum  of  sea  surface  height  variability. 


Area 

Latitude 

Longitude 

R\  (km) 

Cl-max  (Cm/s) 

T\-min  (days) 

1 

11.9-12.1°N 

180-240°  E 

91 

-18.6 

71 

2 

18.9-19. 1°N 

210-250°E 

58 

-7.4 

116 

3 

18.9-19. 1°N 

150-200°E 

67 

-9.7 

100 

4 

21.9-22. 1°N 

145-195°E 

60 

-7.5 

114 

5 

27.9-28. 1°N 

150-210°E 

47 

-4.4 

154 

6 

33. 4-33. 6°  N 

150-190°E 

38 

-2.7 

200 

7 

37.9-38. 1°N 

150-190°E 

30 

-1.6 

270 

Table  8-1:  Areas  chosen  to  analyze  the  frequency/wavenumber  spectra.  Zonal 
and  meridional  ranges,  first  baroclinic  Rossby  radius  of  deformation  Ri  (Chelton  et 
ah,  1998),  theoretical  fastest  phase  speed  of  the  first  baroclinic  long  Rossby  wave 
Ci— max  =  —fiR-i,  and  the  minimum  period  of  the  first  baroclinic  Rossby  waves  in  each 
region  7\_mjn  =  27t /  cf\—max- 
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8.4  Results 


Area  4. 

Consider  Area  4,  northwest  of  the  Hawaiian  chain.  Figure  8-3  depicts  positions  for 
which  altimetric  time  series  were  generated,  thus  defining  the  array.  Fig.8-4a  shows 
the  so-called  beam-pattern  of  the  array  and  which  can  be  interpreted  as  the  apparent 
wavenumber  spectrum  which  would  be  centered  on  wavenumber  (kojo)  if  there  were 
a  monochromatic  wave  of  this  wavenumber  at  a  fixed  frequency  crossing  the  array. 
Notice  that  only  a  very  narrow  meridional  wavenumber  band  is  shown  because  of  the 
lack  of  resolution  in  the  north-south  direction.  Figures  8-4b-f  show  the  result  of  the 
analysis  at  frequencies  plotted  as  logio  [$  (a,  k,  l )  /<3> max  {a->  k,  /)]  with  a  *  denoting  the 
position  of  the  computed  maximum.  As  expected,  the  energy  density  is  generally 
much  larger  for  negative  wavenumbers  (westward  phase  velocity)  than  for  eastward 
going  motions.  There  is  a  finite  background  energy  at  all  wavenumbers  upon  which 
the  peaks  are  superimposed.  This  background  continuum,  which  generally  has  more 
energy  than  is  found  in  the  plane- wave  peak,  requires  separate  study. 

Because  the  analysis  is  essentially  a  two-dimensional  (a,  k )  one,  and  the  1  =  0 
waves  are  the  fastest  ones,  Fig.  8-5  depicts  the  dispersion  relations  for  the  barotropic 
and  first  baroclinic  modes  in  Area  4,  with  1  =  0  (first  baroclinic  Rossby  radii  are 
from  Chelton  et  ah,  1998).  The  observed  spectral  peaks  have  been  transferred  from 
Fig.  8-4  to  Fig.  8-5,  and  an  uncertainty  estimate  is  generated  from  the  confidence 
limits  derived  by  Capon  and  Goodman  (1970).  An  indication  of  the  relative  strength 
of  the  peak  is  provided  by  the  circle  diameter.  The  only  observed  spectral  peaks 
are  on  the  longwave  branch  (wavenumbers  below  the  wavenumber  corresponding  to 
the  maximum  frequency  for  mode  1).  With  one  exception,  at  a  frequency  near  0.008 
cycles/day  (period  123  days),  the  observed  peaks  are  indistinguishable  within  error 
estimates  from  the  conventional  linear  dispersion  curves.  For  the  comparatively  weak 
motions  lying  at  frequencies  a  >  <7i_max  the  motions  cannot  be  explained  as  linear 
baroclinic  waves.  The  fraction  of  the  energy  lying  in  the  observed  peak  varies  from 
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about  43%  near  862  days  period  declining  monotonically  to  about  11%  at  157  days 
period. 

The  zonal  group  velocity  vanishes  at  a  =  on-max  and  under  some  circumstances 
(e.g.,  Wunsch  and  Gill,  1976)  one  anticipates  an  energy  density  maximum  here  in  a 
form  of  resonance.  We  see  no  evidence  for  such  motions. 

Phase  Speed 

By  definition,  the  (zonal)  phase  speed  of  a  wave  is  c  =  a/k,  the  slope  of  the 
dispersion  curve.  For  those  frequencies  and  wavenumbers  in  Fig.  8-5  lying  on  or 
near  the  dispersion  curves,  the  observed  and  theoretical  phase  speeds  are  in  accord 
within  the  uncertainty  limits.  At  higher  frequencies  (or  absolute  wavenumbers),  the 
observed  frequencies  are  systematically  higher  than  the  first  baroclinic  dispersion 
curve  would  permit.  The  average  slope  of  the  empirical  dispersion  curve,  is  higher 
than  the  theoretical  mean  slope  and  one  thus  has  a  positive  phase  speed  bias.  If  we 
rule  out  the  possibility  that  the  waves  are  barotropic  motions  with  too  low  a  phase 
speed,  then  one  could  fairly  conclude  that  these  waves  are  moving  faster  than  theory 
permits. 

Equally  appealing  however,  would  be  their  description  as  motions  not  simply 
consistent  with  linear  dynamics.  Many  sources  of  such  motions  are  known,  in¬ 
cluding  wind-forcing,  instabilities,  non-linear  self-interactions  etc.  On  smaller  scales 
(mesoscale),  the  existence  of  baroclinic  motions  not  satisfying  linear  dispersion  rela¬ 
tions  has  been  known  for  a  long  time  (e.g.,  The  MODE  Group,  1978).  Thus,  if  all  the 
motions  lying  within  a  finite  frequency  band  are  used  to  compute  an  apparent  overall 
phase  velocity,  irrespective  of  the  physics  present,  there  is  clearly  a  bias  toward  large 
values.  But  the  separation  into  distinct  frequency  and  wavenumber  regions  shows 
that  much  of  the  motion  is  not  in  conflict  with  the  simplest  linear  theory. 
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frequency  (cycles  per  day) 


Figure  8-5:  Frequency-zonal  wave  number  spectra  of  region  4.  The  peak  in  each 
frequency  band  in  Fig.8-4  is  indicated  by  a  dot  here.  The  magnitude  of  each  peak 
is  denoted  by  the  size  of  the  circle  around  that  peak.  The  dash-dot  and  solid  curves 
denote  the  barotropic  and  first  baroclinic  Rossby  wave  dispersion  curves  with  l  =  0, 
respectively.  At  the  lowest  frequencies,  the  barotropic  curve  is  visually  indistinguish¬ 
able  from  k  =  0.  Uncertainty  estimates  are  derived  from  the  95%  confidence  interval 
(Capon  and  Goodman,  1970)  for  the  peak  values.  The  zonal  wavenumber  interval  rep¬ 
resents  those  wavenumbers  where  the  background  energy  surrounding  the  apparent 
peak  is  larger  than  or  equal  to  the  lower  limit  of  the  confidence  interval. 
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Other  Areas. 

The  frequency/zonal-wavenumber  spectra  for  the  remaining  areas  analyzed  here 
are  shown  in  Fig.  8-6.  In  Area  1  (Fig.  8-6a)  ,  most  peaks  are  located  at  k  —  0.  These 
motions  have  zonal  scales  larger  than  the  array  aperture  and  can  be  due  to  large-scale 
barotropic  motions  or  steric  forcing  on  large  zonal  scales  (annual  cycle  heating)  or 
both.  Only  two  peaks  are  not  located  at  k  =  0  and  they  agree  very  well  with  the 
dispersion  curve  for  first  baroclinic  Rossby  waves,  but  carry  only  about  10%  of  the 
energy.  Here  again,  a  lumped  apparent  phase  speed  would  be  biassed  high. 

As  shown  in  Fig.  8-1,  the  latitudes  of  Areas  2  and  3  are  the  same,  east  and  west 
respectively  of  the  Hawaiian  chain.  The  frequency  zonal-wavenumber  spectra  for  these 
areas  are  shown  in  Fig.  8-6b,c.  As  indicated  in  Table  8-1,  long  first  baroclinic  mode 
waves  propagate  faster  in  the  western  parts  of  the  basin.  This  zonal  change  of  phase 
speed  is  a  result  of  the  deepening  of  the  thermocline  in  the  west  (Chelton  et  ah,  1998). 
At  the  annual  frequency  band,  Area  3  is  dominated  by  a  first  baroclinic  Rossby  wave, 
whereas  Area  2  is  dominated  by  a  signal  indistinguishable  from  zonal  uniformity 
( k  =  0),  and  which  is  presumably  the  steric  contribution  alluded  to  above.  One 
explanation  for  this  difference  (e.g.,  Chelton  and  Schlax,  1996)  is  that  the  Hawaiian 
Ridge  generates  a  westward  moving  train  of  annual  period  Rossby  waves.  Fig.  8-6b 
and  8-6c  show  that  at  low  frequencies  the  observed  spectral  peaks  are  located  either 
very  close  to  the  baroclinic  dispersion  curve  (around  20%  of  the  total  energy)  or  at 
the  axis,  k  —  0.  At  frequencies  higher  than  <7i_mai,  the  observed  spectral  peaks  are 
situated  between  the  barotropic  and  first  baroclinic  Rossby  wave  dispersion  curves. 

As  the  latitude  increases,  the  maximum  frequency  and  the  phase  speed  of  the 
lowest  baroclinic  mode  (the  fastest  baroclinic  mode)  decreases.  Whatever  baroclinic 
motion  is  present  above  <7i_max,  it  cannot  be  a  linear  flat-bottom  mode.  Even  where 
the  waves  are  present  in  linear  theory,  they  become  much  more  susceptible  to  the 
various  effects  described  by  Killworth  et  al.  (1997)  and  others. 

White  (1977)  inferred  the  existence  of  annual-period  baroclinic  Rossby  waves  in 
the  tropical  Pacific  and  he  suggested  that  they  are  generated  by  the  annual-period 
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wind  forcing.  Here  we  find  that  only  Areas  3  and  4  show  such  waves  significantly 
present.  Other  regions  are  dominated  by  energy  indistinguishable  from  k  =  0  and 
which  may  be  spectral  leakage  from  the  annual  period. 

Within  the  equatorial  band  from  10°N  to  10°S,  CS  (their  Fig.  5)  found  that  appar¬ 
ent  westward  phase  speeds  were  smaller  than  the  theoretical  value.  At  these  latitudes, 
one  expects  that  higher  modes  should  become  more  prominent,  and  this  is  observed 
in  current  meter  records  (e.g.,  Wunsch,  1997)  and  in  data  from  profiling  devices  (e.g., 
Eriksen,  1981).  The  unaccounted-for-presence  of  high  modes  is  consistent  with  a  neg¬ 
atively  biassed  value  for  C\.  Furthermore,  around  the  equator  eastward  propagating 
Kelvin  and  mixed-Rossby-gravity  wave  modes  become  possible.  The  equatorial  band 
is  thus  worthy  of  special  attention,  and  it  will  not  be  further  discussed  here. 

The  Short- Wavelength  Limit 

High-resolution  methods  permit  one  to  use  irregularly  spaced  arrays  to  seek  waves 
whose  wavelengths  are  set  by  the  Nyquist  wavelength  of  the  shortest  lateral  separa¬ 
tion.  Robustness  of  the  results  is  dependent  upon  signal-to-noise  ratios  and  the  actual 
array  configuration.  An  attempt  was  made  in  Area  4  to  determine  the  wavenum¬ 
ber  spectrum  at  wavenumbers  beyond  those  corresponding  to  a\-max,  the  short- 
wavelength  branch.  No  statistically  significant  peaks  above  the  background  energy 
could  be  found,  and  we  will  not  discuss  these  results  any  further. 

8.5  Discussion 

We  find  a  significant  fraction  of  the  energy  present  in  the  North  Pacific  Ocean  is 
consistent  with  the  linear  Rossby  wave  dispersion  relation  Eq.  (2.36).  The  remaining 
energy  is  not  so  consistent,  and  believe  that  the  results  of  CS  are  readily  explained 
through  the  general  mechanisms,  if  not  all  the  details,  described  by  Qiu  et  al.  (1997), 
Killworth  et  al.  (1997)  and  the  wider  literature  on  generalizations  to  more  complex 
conditions  of  the  linear  theory.  One  has  in  part  an  interference  pattern  of  free  motions 
consistent  with  linear  physics,  superimposed  upon  forced  motions  of  various  types  . 
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In  any  fixed  wavenumber  band,  there  exists  some  motions  with  frequencies  lying 
above  the  dispersion  curve  for  first  baroclinic  modes.  These  motions  produce  lumped 
westward  phase  speeds  which  are  generally  too  high  for  the  linear  theory.  It  is  almost 
a  semantic  distinction  as  to  whether  the  motion  is  best  described  as  linear  Rossby 
waves  with  the  “wrong”  frequency,  or  simply  as  oceanic  motions  not  satisfying  the 
dispersion  relation.  At  most  latitudes,  there  is  significant  energy  which  does  appear 
consistent  within  present  error  estimates  with  the  linear  dispersion  relationship. 

The  question  raised  by  CS  as  to  why  c„  is  too  large,  now  becomes  instead  the 
question  of  why  frequencies  become  preferentially  too  high  as  wavenumber  magnitudes 
increase?  Presumably  the  answer  to  this  question  is  some  combination  of  the  answer 
given  by  Killworth  et  al.  (1997),  the  mode  coupling  of  the  barotropic  and  baroclinic 
modes  known  to  occur  in  some  regions  (e.g.,  Wunsch,  1997),  and  the  space/time 
structure  of  the  atmospheric  forcing.  Many  reasons  can  be  found  for  the  existence 
of  motions  inconsistent  with  the  first  baroclinic  mode  of  a  resting  flat-bottom  ocean, 
including  the  presence  of  the  barotropic  mode,  higher  baroclinic  modes,  advection, 
topography  and  forcing. 

At  higher  latitudes,  as  the  waves  slow,  and  ultimately  disappear  at  frequencies 
accessible  with  existing  data  sets,  non-linear  effects  as  well  as  background  distortions 
of  the  potential  vorticity  gradient  (e.g.  Killworth,  et  al.,  1997)  can  dominate  the 
spectrum  and  eventually  there  is  no  evidence  for  linear  waves  at  all. 
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Chapter  9 


Conclusions  and  discussion 


This  thesis  work  can  be  summarized  as  producing  an  approximate  basic  spectral  de¬ 
scription  of  low  frequency  oceanic  variability  by  using  a  simple  dynamical  model  with 
various  observations.  First,  this  study  shows  that  a  single  function  of  &(k,  l,  uj,  n,  A,  (f>) 
does  exist,  and  all  measurements  can  be  simplified  as  a  certain  function  of  $(&,  l,  uj,  n ,  A,  (/>). 
It  implies  that  various  data  are  consistent  with  linear  equations  of  fluid  dynamics  to 
first  order  approximation.  This  work  provides  estimates  of  how  the  energy  of  low 
frequency  oceanic  variability  depends  on  two  components  of  horizontal  wavenumbers, 
frequency,  vertical  mode  and  geography,  and  what  percentage  of  energy  propagates 
westward  and  eastward.  The  only  geography-dependent  part  is  the  energy  level.  It 
also  provides  estimates  of  some  physical  variables  in  which  no  observations  are  avail¬ 
able,  such  as  the  wavenumber  spectrum  of  horizontal  velocities.  With  the  aid  of 
Gauss-Markov  estimation  method,  the  spectrum  is  applied  to  the  design  of  observa¬ 
tional  strategy.  If  the  regional  model  $(k,  l,  uj,  n,  A,  <f>)  is  normalized  by  the  geography- 
dependent  factor  /(</>,  A),  we  will  obtain  a  universal  spectral  model  ${k,l,uo,n).  The 
vertical  structure  is  governed  by  the  barotropic  and  first  two  baroclinic  modes,  and 
the  frequency  and  wavenumber  structure  of  each  mode  is  displayed  in  figure  9-1. 

The  low  frequency  variability  is  a  broad-band  process  and  Rossby  waves  are  par¬ 
ticular  parts  of  it.  The  sea  surface  height  measurements  in  the  North  Pacific  are 
analyzed,  and  it  is  found  that  about  10-40%  of  the  total  variance  in  each  frequency 
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band  are  indistinguishable  from  the  simplest  theoretical  Rossby  wave  description. 
Toward  higher  latitudes,  both  amplitudes  and  phase  speeds  of  linear  waves  decrease. 
Non-equatorial  latitudes  display  some  energy  with  frequencies  too  high  for  consis¬ 
tency  with  linear  theory;  this  energy  produces  a  positive  bias  if  a  lumped  average 
westward  phase  speed  is  computed  for  all  the  motions  present. 

The  model  produced  here  is  not  final  and  will  become  obsolete  with  the  coming 
of  new  observations.  For  example,  the  vertical  structure  of  kinetic  energy  has  been 
obtained  by  Wunsch  (1997)  from  current  meter  data;  however,  each  mode  cannot  be 
completely  represented  by  the  energy  level  alone.  Each  vertical  mode  is  a  function  of 
two  components  of  horizontal  wavenumber  and  frequency.  So  far,  complete  frequency 
and  wavenumber  structure  of  barotropic  and  baroclinic  modes  of  any  field  still  does 
not  exist  except  that  a  few  current  meter  moorings  can  provide  an  estimate  of  it;  how¬ 
ever,  the  resolution  and  reliability  of  such  estimate  are  very  low.  The  coming  Gravity 
Recovery  and  Climate  Experiment  (GRACE)  project  will  enable  oceanographers  to 
measure  ocean  bottom  pressure  fluctuations  from  space  with  extraordinary  accuracy 
through  measuring  the  time-dependent  gravity.  Simultaneous  measurements  of  sea 
surface  height  and  ocean  bottom  pressure  will  greatly  enhance  our  knowledge  of  the 
ocean  circulation.  Among  other  things,  this  can  help  us  to  estimate  how  ocean  vari¬ 
ability  is  partitioned  between  barotropic  and  baroclinic  motions  and  the  space  and 
time  scales  of  barotropic  motions,  which  could  be  used  to  improve  the  energy  density 
model  produced  in  this  thesis.  The  present  model  has  the  following  assumptions: 

(1)  observations  at  different  times  can  be  reasonably  related,  which  implies  the 
oceanic  variability  is  stationary  in  time, 

(2)  the  process  is  linear, 

(3)  the  basic  state  is  at  rest, 

(4)  the  bottom  of  the  ocean  is  flat, 

(5)  there  is  no  forcing,  and 

(6)  the  buoyancy  frequency  N(z )  is  universal  and  can  be  modeled  by  an  idealized 
exponential  form. 
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The  empirical  model  produced  in  this  thesis  can  only  be  applied  to  the  mid-ocean 
away  from  the  equator  and  boundaries.  Regions  close  to  the  equator  and  western 
boundary  have  distinct  dynamics.  Near  the  western  boundary  currents,  the  barotropic 
mode  dominates  (Wunsch,  1997);  the  effect  of  mean  flow  cannot  be  neglected  and 
there  exist  significant  differences  between  the  energy  level  of  zonal  and  meridional 
components  of  horizontal  velocity.  The  frequency  and  wavenumber  spectral  shape 
near  the  western  boundary  is  different  from  that  in  the  interior  ocean  (Stammer, 
1997).  In  the  areas  close  to  the  equator,  high  modes  become  more  important  (Wunsch, 
1997;  Eriksen,  1981),  and  the  geostrophic  balance  doesn’t  hold.  As  discussed  in 
chapter  2,  my  model  cannot  represent  temperature  variations  in  the  mixed  layer  and 
at  sea  surface.  Most  data  used  in  this  study  span  a  few  years.  The  frequency  and 
wavenumber  spectrum  for  very  low  frequencies — decadal  and  even  longer  period — is 
unclear.  The  seasonal  cycle  requires  further  study  also.  As  displayed  in  chapter  5, 
there  is  a  striking  peak  associated  with  the  seasonal  cycle  of  the  frequency  spectrum 
of  sea  surface  height.  However,  the  seasonal  cycle  is  not  significant  on  the  frequency 
spectrum  of  horizontal  velocity  and  temperature  below  the  mixed  layer. 

Why  <b(k,l,  u,  n)  has  a  particular  universal  spectral  shape  is  of  great  interest 
from  the  theoretical  point  of  view.  The  generation  mechanism  of  ocean  variability 
includes  atmospheric  wind  and  thermal  forcing,  barotropic  and  baroclinic  instability, 
and  the  interaction  of  mean  flow  with  topography.  The  lowest-mode  character  of  low 
frequency  variability  is  attributed  to  quasi-geostrophic  nonlinear  interactions,  which 
drives  the  motion  toward  larger  scales  both  in  the  vertical  and  horizontal  direction 
(Charney,  1971;  Rhines,  1975;  Fu  and  Flierl,  1979).  In  the  ocean  and  atmosphere, 
wavenumber  spectral  shape  is  often  rationalized  by  geostrophic  turbulence  theory. 
In  the  absence  of  forcing  and  dissipation,  two-dimensional  flow  is  peculiar  in  that 
total  kinetic  energy  and  enstrophy  are  conserved.  Therefore,  kinetic  energy  must 
be  transferred  simultaneously  from  small  to  larger  scales  and  large  to  smaller  scales 
through  nonlinear  interactions  between  different  scales  of  motion  (Fjortoft  1953). 
This  led  Krainchnan  (1967)  to  postulate  that  there  exists  two  distinct  “inertial” 
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subranges  of  two-dimensional  turbulence.  The  large-scale  subrange  has  a  kinetic 
energy  spectrum  of  k~ 3/5,  while  the  kinetic  energy  spectrum  of  small-scale  subrange 
is  proportional  to  Ar3.  Charney  (1971)  also  obtained  a  —3  law  at  the  tail  of  the  kinetic 
energy  spectrum  from  a  quasi-geostrophic  model  and  conjectured  that  in  this  region 
there  would  be  equipatition  between  the  two  components  of  the  kinetic  energy  and 
the  potential  energy.  Rhines  (1979)  further  argued  that  kinetic  energy  should  exhibit 
a  maximum  at  the  wavelength  close  to  the  “Rhines  scale” ,  l  oc  ^Ju/ /?,  where  u  is  the 
root-mean-square  of  the  fluctuating  speed.  The  results  from  numerical  simulation 
are  quite  similar  (Lilly,  1969;  McWilliams  and  Chow,  1981).  For  the  ocean,  two- 
dimensional  turbulence  theory  only  offers  an  explanation  for  the  spectral  shape  at 
wavelengths  smaller  than  400  km.  It  is  unclear  why  the  wavenumber  spectra  of  sea 
surface  height  and  temperature  are  proportional  to  A;-1/2  at  wavelengths  longer  than 
400  km.  The  frequency  spectral  shape  is  explained  as  an  integral  response  of  the 
ocean  to  continuous  random  forcing  of  the  atmosphere  (Hasselman,  1976).  Using 
a  conceptual  model  with  atmospheric  forcing,  Hasselman  (1976)  and  Frankignoul 
and  Hasselman  (1977)  found  that  although  the  frequency  spectrum  of  atmospheric 
forcing  is  white  at  periods  longer  than  one  week,  the  frequency  spectrum  of  the  ocean 
should  be  red  and  proportional  to  cu-2  at  periods  from  a  few  days  to  about  150  days 
because  the  ocean  acts  as  a  first  order  Markov  integrator  of  atmospheric  input.  At 
periods  longer  than  a  year,  the  frequency  spectrum  of  the  ocean  turns  white,  because 
negative  feedbacks  become  important.  Frankignoul  and  Muller  (1979)  and  Willebrand 
et  al.  (1980)  further  investigated  the  response  of  the  ocean  to  stochastic  forcing  by 
the  atmosphere  using  linear  quasi-geostrophic  model.  Willebrand  et  al.  (1980)  found 
that  the  ocean  response  is  essentially  barotropic  at  periods  between  the  inertial  period 
and  about  300  days.  Most  of  the  theories  only  provide  qualitative  explanation  with 
focus  on  predicting  the  wavenumber  and  frequency  spectral  shape.  Even  the  state- 
of-art  general  circulation  model  cannot  predict  both  spectral  shape  and  energy  levels 
correctly  (Stammer  et  al.,  1996).  More  work  needs  to  be  done  along  these  lines  in 
the  future. 
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Appendix  A 


Rossby  waves  and  the  mid- Atlantic 
ridge 

A.l  Introduction 

Rossby  waves  are  large  scale  and  low  frequency  phenomena  in  the  ocean.  The  de¬ 
tection  of  Rossby  waves  is  only  possible  with  a  global  data  set.  Satellite  altimetry 
enables  oceanographers  to  monitor  the  global  ocean  variability  continuously  for  years 
for  the  first  time.  The  Rossby  waves  have  been  identified  by  many  oceanographers 
from  the  satellite  altimeter  data  (e.g.,  Chelton  and  Schlax,  1996;  Wang  et  al,  1998; 
Zang  and  Wunsch,  1999).  The  generation  mechanism  of  Rossby  wave  includes:  (1) 
atmospheric  wind  forcing  and  buoyancy  forcing  acting  on  the  sea  surface,  (2)  either 
along-shore  wind  stress  or  coastally  trapped  Kelvin  waves  along  the  eastern  bound¬ 
ary,  and  (3)  the  interaction  between  the  mean  flow  and  topography.  The  generation 
of  Rossby  waves  by  wind  stress  has  received  much  attention  in  the  literature  (e.g., 
Veronis  and  Stommel,  1956;  Anderson  and  Gill,  1975;  White,  1977;  and  Qiu  et  al., 
1997).  Chelton  and  Schlax  (1996)  found  that  only  in  low  latitudes,  could  boundary- 
gen  erated  Rossby  waves  propagate  across  the  ocean  basin.  Towards  high  latitudes, 
the  boundary-originated  waves  were  trapped  near  the  east  coast.  This  feature  was 
also  found  in  numerical  ocean  modeling  (Qiu  et  al.,  1997).  Recent  theoretical  studies 
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suggest  that  the  ocean  general  circulation  is  strongly  influenced  by  mid-ocean  ridge 
system,  which  acts  as  an  incomplete  barrier  to  the  ocean  circulation  (Pedlosky  et  ah, 
1997;  Pedlosky  and  Spall,  1999).  One  open  question  is  whether  the  Rossby  waves 
generated  at  the  eastern  boundary  can  cross  the  mid-ocean  ridge.  If  the  anoma¬ 
lies  originated  from  the  eastern  boundary  can  be  transmitted  across  the  mid-ocean 
ridge  as  westward-going  Rossby  waves,  they  will  affect  the  western  boundary  currents; 
thus  influence  the  climate  in  other  regions.  Using  a  quasi-geostrophic  model,  Barnier 
(1988)  investigated  the  influence  of  the  mid- Atlantic  ridge  on  first  baroclinic  Rossby 
waves  generated  by  seasonal  wind  fluctuations  at  the  eastern  boundary.  He  found 
that  the  ridge  blocked  the  waves  coming  from  the  eastern  boundary  but  generated 
new  waves  whose  phase  vector  was  almost  normal  to  the  ridge  crest  west  of  the  ridge. 

In  general,  the  propagation  of  Rossby  waves  can  be  represented  by  three  different 
ways.  The  direct  method  is  to  present  the  measurements  in  space  and  time  domain: 
the  so-called  hovmoller  diagram.  Strictly  speaking,  this  method  is  useful  only  when 
the  signal  is  made  up  of  a  single  wave.  If  the  data  is  composed  of  waves  with  dif¬ 
ferent  frequencies  and  wavenumbers,  the  spectral  method  is  required.  The  spectral 
method  enables  one  to  quantitatively  determine  the  wavenumbers  and  frequencies  of 
dominant  signals  in  the  data.  It  also  makes  it  easy  to  compare  observations  with 
theory.  However,  the  spectrum  provides  no  information  about  time  and  space.  If 
the  signal  undergoes  abrupt  changes  over  time  or  space,  the  Fourier  method  is  not 
suitable  to  detect  these  abrupt  changes.  The  third  way  of  presenting  propagation 
of  Rossby  waves  is  wavelet  transforms.  The  main  advantage  of  wavelets  is  that  they 
contain  information  about  both  frequency  (wavenumber)  and  time  (space).  A  wavelet 
transform  is  a  powerful  tool  for  the  study  of  a  rapidly  changing  wave  field  during  wave 
generation,  growth  or  decay. 

In  this  chapter,  the  array  signal  processing  technique  will  be  used  with  wavelet 
method  to  study  how  the  properties  of  Rossby  waves  change  across  the  mid- Atlantic 
ridge.  First,  the  array  signal  processing  technique  will  be  applied  to  TOPEX/POSEIDON 
satellite  altimeter  data  to  find  out  the  period  and  wavelength  of  dominant  Rossby 
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waves.  Then,  the  wavelet  technique  will  be  used  to  investigate  how  the  property  of 
dominant  Rossby  waves  changes  with  space  and  time.  The  objective  here  is  to  find 
out  whether  Rossby  waves  west  of  major  ocean  ridges  are  generated  locally  by  the 
topography  or  originate  from  the  eastern  boundary  and  amplify  over  the  mid-ocean 
ridge. 

A. 2  Data  and  data  processing 

TOPEX/POSEIDON  altimeter  along  track  data  from  cycle  11  to  212  (December  31, 
1992  to  June  15,  1998)  is  used  here.  In  order  to  study  the  effect  of  the  mid-ocean 
ridge  on  Rossby  waves,  one  particular  box  in  the  North  Atlantic,  which  crosses  the 
mid-Atlantic  ridge,  is  chosen  to  study.  The  spatial  distribution  of  measurements 
within  the  box  is  fairly  regular  (Figure  A-l).  TOPEX/POSEIDON  provides  global 
sea  surface  height  measurements  every  9.91  days.  The  measurements  within  one  orbit 
cycle  are  assumed  to  be  synchronous.  First  the  time  mean  of  the  sea  surface  height 
at  each  spatial  point  is  removed,  then  a  three  block  non-overlapping  time  average  was 
done  to  remove  motions  with  periods  shorter  than  59.5  days. 

A. 3  Wavelet  analysis 

Wavelet  method  was  developed  in  last  decade,  since  then  it  has  been  applied  to 
the  study  of  various  fields  such  as  decadal  variability,  atmospheric  turbulence  and 
ocean  wind  waves  (see  Kumar  and  Foufoula-Georgiou,  1997  for  a  review).  Wavelet 
method  provides  quantitative  information  about  not  only  the  wavenumber,  frequency 
and  phase  speed  of  propagating  waves  but  also  spatial  and  temporal  modulation  of 
energy  amplitude.  Wavelet  transform  and  the  relationship  between  wavelet  transform 
and  Fourier  transform  are  briefly  summarized  here. 
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A. 3.1  Wavelet  transform 


The  continuous  wavelet  transform  S(b,a )  of  a  real  signal  f(t)  with  respect  to  the 
mother  wavelet,  g(t),  is  defined  as 
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represents  a  set  of  functions,  that  is,  wavelets,  which  are  constructed  from  the  mother 

wavelet  g(t)  by  translating  and  dilating.  Here  a  is  a  scale  parameter,  b  is  a  local 

parameter  and  *  stands  for  complex  conjugate.  The  normalizing  constant  -K*  is 
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chosen  so  that  gab(t)  has  the  same  energy  for  all  scales  a: 
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The  inverse  wavelet  transform  is  given  by 
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where  G(u>)  is  the  Fourier  transform  of  g(t), 
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Using  Parseval’s  theorem,  one  can  write  the  function  S(b,  a)  in  terms  of  F(u)  and 
G(u>),  which  are  the  Fourier  transforms  of  f(t )  and  g(t ): 


S(b,a)  =  \/a  f  F(u))G*(au))e^u duj. 
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(A. 7) 


164 


A. 3. 2  Wavelet  and  Fourier  transform 


The  advantages  of  wavelet  transform  have  been  discussed  in  the  introduction;  how¬ 
ever,  wavelet  analysis  has  its  own  limitations.  Fourier  transform  is  orthogonal  and 
is  a  standard  procedure  to  solve  partial  differential  equations  analytically.  Moreover, 
because  most  theory  on  waves  is  expressed  in  terms  of  frequency  and  wavenumber, 
it  is  straightforward  to  compare  spectrum  of  the  observations  with  theory.  However, 
continuous  wavelet  transform  is  not  orthogonal;  thus  it  cannot  be  used  to  solve  partial 
differential  equations  analytically.  In  addition,  because  the  wavelet  transform  is  not 
orthogonal,  the  wavelet  transform  may  produce  some  false  correlations  (see  Farge, 
1992).  One  must  beware  of  this  redundancy  when  applying  the  wavelet  transform  to 
the  actual  data. 

As  indicated  by  equation  (A. 7),  the  wavelet  transform  depends  on  both  the  signal 
and  the  choice  of  wavelet.  The  choice  of  wavelet  is  contingent  on  what  kind  of 
information  one  wants  to  extract  from  the  data.  The  Morlet  wavelet  has  been  widely 
used  in  atmosphere  and  oceanography  because  Morlet  wavelet  is  suitable  for  the 
study  of  wave  propagation  and  it  is  easy  to  compare  the  Morlet  wavelet  transform 
with  Fourier  transform. 

The  Morlet  wavelet  is  defined  as 

git)  =  t T^e-^e-*2/2,  (A.8) 

which  is  a  periodic  function  modulated  by  a  Gaussian  envelop.  It  is  complexed- valued; 
therefore,  it  enables  one  to  extract  the  information  about  the  amplitude  and  phase  of 
the  data  being  analyzed.  In  order  to  satisfy  the  admissibility  condition,  it  is  necessary 
to  take  c  >  5  (Farge,  1992). 

In  order  to  compare  wavelet  transform  of  observations  with  wave  theory,  one  needs 
to  understand  the  relationship  between  the  wavelet  transform  and  Fourier  transform, 
particularly,  what  is  the  relationship  between  wavelet  scales  and  Fourier  periods?  and 
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what  is  the  relationship  between  wavelet  amplitude  and  Fourier  amplitude?  Meyers 
et  al.  (1993)  demonstrated  that  in  the  case  of  Morlet  wavelet,  the  relationship  be¬ 
tween  Wavelet  scales  and  Fourier  periods  could  be  derived  analytically.  They  showed 
that  the  scale,  oo,  where  the  Morlet  wavelet  transform  of  a  periodic  signal  Aoel2nt/T° 
reached  its  maximum  was  related  to  the  period  of  the  signal,  T0  through 


c+p  +  c2)1/2 
a°  —  l  7Z - Fo, 


47T 


(A.9) 


where  c  is  the  parameter  of  Morlet  wavelets.  If  one  chooses 


c  —  2n  — 
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then  the  scale  of  wavelet  transform  equals  the  period  of  the  signal: 


(A. 10) 


oo  =  r0.  (A.ii) 

It  is  ready  to  demonstrate  that  the  corresponding  relationship  between  the  maximum 
of  the  square  amplitudes  of  the  wavelet  transform  S(b,  a)  and  the  amplitude  of  the 
wave  A0  is 

a)lmax  =  27r“1/2e1/(47r)Vo|2a.  (A.12) 

Therefore,  the  amplitude  of  wavelet  transform  not  only  depends  on  the  amplitude  of 
the  wave  but  also  depends  on  the  scale  (period)  of  the  wave.  In  order  to  compare  the 
true  energy  at  different  scales  unbiasly,  it  is  desirable  to  present  \S(b,a)\2/a,  rather 
than  |5(f>,  a)|2,  which  is  widely  used  in  the  literature. 
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A.4  Application 

A. 4.1  Frequency  and  wavenumber  spectrum 

The  three-dimensional  spectrum  of  sea  surface  height  anomalies  is  obtained  by  using 
the  conventional  array  signal  processing  method.  Non-overlap  averaging  over  three 
adjacent  frequency  bands  is  performed  so  that  statistical  stable  estimates  are  ob¬ 
tained.  The  three-dimensional  spectrum  is  plotted  in  figure  A-2  as  logw  [<3>  (a;,  k ,  l )  / <J>m 
with  a  *  denoting  the  position  of  the  computed  maximum.  As  displayed  in  figure  A-2, 
the  peak  at  each  frequency  band  occurs  either  at  the  axis  of  negative  zonal  wavenum¬ 
bers  or  at  the  origin.  There  is  more  energy  going  westward  just  as  expected.  The 
annual  frequency  band  is  dominated  by  stationary  waves,  which  is  presumably  due  to 
uniform  atmospheric  forcing.  At  the  periods  of  178  and  245  days,  the  fraction  of  the 
energy  lying  in  the  observed  peak  is  about  15%,  which  is  higher  than  that  at  other 
periods.  The  zonal  wavelength  of  the  observed  peak  at  178  and  245  days  is  700  and 
930  km,  respectively.  The  observed  phase  speeds  at  these  two  periods  are  very  close, 
roughly  4.4  cm/s. 

The  averaged  variance  of  sea  surface  height  anomaly  in  the  area  studied  here  is 
about  45  cm2,  among  which  about  7  cm2  is  contained  in  eastward-going  motions, 
22  cm2  in  westward- going  motions  and  16  cm2  in  stationary  motions.  Figure  A- 
3a  shows  the  frequency  spectra  of  eastward-going,  westward-going  and  stationary 
motions.  There  is  a  pronounced  peak  at  the  period  of  one  year  in  the  frequency 
spectrum  of  stationary  motions.  The  energy  of  westward-going  motions  is  higher  than 
the  energy  of  eastward-going  motions  at  almost  all  frequency  bands.  The  frequency 
spectrum  of  eastward-going  motions  is  red,  and  the  frequency  spectrum  of  westward¬ 
going  motions  is  dominated  by  motions  with  periods  from  150  days  to  300  days. 
The  frequency  spectrum  of  westward-going  motions  is  red  at  periods  shorter  than 
150  days  and  is  blue  at  periods  longer  than  300  days.  The  wavenumber  spectra 
of  eastward-going  and  westward-going  motions  are  displayed  in  figure  A-3b.  Most 
striking  property  in  figure  A-3b  is  that  there  is  a  peak  at  the  wavelength  about  750 
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km  in  the  wavenumber  spectrum  of  westward-going  motions. 

Because  the  meridional  range  of  the  area  in  figure  A-l  is  so  small  that  the  analysis 
is  a  two-dimensional  one  in  essence.  The  barotropic  and  baroclinic  dispersion  curves 
for  classic  Rossby  theory  with  l  =  0  are  displayed  in  figure  A-4.  Each  dot  in  figure 
A-4  corresponds  to  the  observed  spectral  peak  at  each  frequency  band  in  figure  A- 
3.  The  size  of  the  circle  stands  for  the  relative  amplitude  of  the  peak.  The  error 
bar  is  produced  according  to  Capon  and  Goodman  (1971).  The  result  in  figure 
A-4  is  consistent  with  those  found  by  Zang  and  Wunsch  (1999)  in  North  Pacific. 
The  observed  spectral  peaks  are  all  located  on  the  longwave  branch  and  the  peaks 
fall  between  the  barotropic  and  baroclinic  dispersion  curves.  The  frequencies  of  the 
observed  peaks  are  systematically  higher  than  those  of  the  first  baroclinic  dispersion, 
and  lower  than  those  of  the  barotropic  dispersion.  Therefore,  the  observed  phase 
speeds  are  between  first  baroclinic  Rossby  wave  speeds  and  barotropic  Rossby  wave 
speeds. 

A. 4. 2  Wavelet  transform 

As  shown  in  figure  A-l,  data  inside  the  box  are  nearly  uniformly-distributed.  Here 
they  are  mapped  into  one  and  half  degree  regular  grid  along  28.3°N  using  linear  inter¬ 
polation  technique.  Since  the  signal  of  Rossby  waves  within  the  periods  between  150 
and  300  days  is  stronger  than  those  at  other  periods,  my  concern  here  is  focused  on 
the  frequency  band  between  1/300  and  and  1/150  cycles  per  day.  The  steric  height 
is  extracted  by  fitting  the  data  with  annual  harmonics  and  is  then  substracted  from 
the  data.  Then,  a  zero-phase  bandpass  filter  is  applied  to  extract  the  signals  with 
periods  between  150  and  300  days  from  the  residuals.  The  Hovmoller  diagram  of  the 
filtered  data  is  shown  in  figure  A-5.  As  displayed  there,  there  is  a  systematic  west¬ 
ward  propagation.  West  of  315°E,  the  amplitudes  of  westward  propagating  signals 
are  much  larger.  Moreover,  the  slope  of  constant  phase  line  is  steeper  east  of  295°E; 
thus,  the  phase  speed  is  higher  in  the  region  close  to  western  boundary.  The  change 
of  phase  speed  with  longitude  results  from  the  fact  that  thermocline  is  shallower 
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in  the  east.  Figure  A-5  also  shows  that  westward  propagation  changes  with  time. 
There  is  an  increase  of  energy  around  900  days.  The  hovemoller  diagram  provides 
only  qualitative  description  of  spatial  and  temporal  variations  wave  propagation.  To 
be  quantitative,  one  can  resort  to  windowed  Fourier  transform  or  wavelet  technique. 
Compared  with  windowed  Fourier  transform,  wavelet  method  is  more  flexible  and  con¬ 
venient.  I  perform  one-dimensional  wavelet  transform  of  the  filtered  sea  surface  height 
data  r](x,  t)  with  respect  to  x.  The  mother  wavelet  used  here  is  Morlet  wavelet  with 
c  =  2tt  —  1/27T.  The  final  product  is  a  three-dimensional  one:  \S(a,x,t)\2/a,  where  a 
is  the  characteristic  scale  (wavelength).  The  temporal  average  of  \S(a,x,t)\2/a  is  dis¬ 
played  in  figure  A-6.  As  displayed  there,  the  variability  has  a  dominant  characteristic 
zonal  scale,  about  730  km  and  is  much  stronger  west  of  the  mid-Atlantic  ridge.  To 
highlight  the  variations  of  energy  amplitude  with  longitude,  the  temporal  average  of 
| S(a,  x,  t)\2/a  with  the  characteristic  scale  a  —  730  km  is  plotted  in  figure  A-7,  which 
shows  that  variability  with  wavelength  about  750  km  increases  rapidly  from  340°E  to 
300°E  and  reaches  its  maximum  near  300°E.  Figure  A-8  displays  the  time  evolution  of 
| S(a,  x,  t)\2/a  with  scale  of  730  km.  There  are  two  obvious  sudden  increases  of  energy 
in  time:  one  near  900  days  and  the  other  near  1600  days.  The  reason  for  these  two 
sudden  changes  requires  further  study,  presumably  associated  with  variations  of  wind 
forcing.  The  corresponding  phase  change  with  time  and  space  is  shown  in  figure  A-9, 
which  indicates  that  variability  propagates  westward  systematically.  The  period  is 
about  200  days  and  the  phase  speed  is  4.2  cm/s.  Figure  A-9  shows  that  although  the 
amplitude  of  the  Rossby  waves  is  much  stronger  west  of  topography,  the  propagation 
of  Rossby  wave  west  of  mid- Atlantic  ridge  originates  from  the  eastern  boundary.  This 
suggests  that  the  Rossby  waves  west  of  mid- Atlantic  ridge  is  generated  along  the  east 
boundary  and  suddenly  amplifies  over  topography. 
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Figure  A-l:  Array  geometry  in  the  area  chosen  to  study.  The  position  of  each  “sen¬ 
sor”  is  marked  by  a  star  and  corresponds  to  a  position  on  the  TOPEX/POSEIDON 


sub-satellite  track. 


1 


0.04 

a> 

-O 

|  0.02 

a» 

> 

co 

I  0 

3 

o 

CO 

J.  -0.02 
o 

-0.04 

-2 

004 

a3 

_Q 

§  0  02 

CD 

> 

CD 

I  0 

3 

o 

CO 

'  -0.02 


beam  pattern 

a 

-0.04 


_ : _ 

61c 

lays 

1 

If 

:  b 

-4  —2 


Figure  A-2a  The  beam  pattern  for  the  array  in  Fig.  1.  Fig.  A-2b-j.,  The  frequency- 
wavenumber  spectra  $(o»,  k,l)  at  periods  T  =  27r/ <x  (in  days  )  (b)  61.  (c)  68.  (d)  75. 
(e)  85.  (f)  98.  (g)  115.  (h)  140.  (i)  178.  (j)  245.  (k)  392.  and  (i)  981  days.  The 
contour  levels  are  the  same  in  each  panel  and  vary  from  0  to  -0.9  in  steps  of -0.3.  The 
units  of  east- west  wavenumber  and  north-south  wavenumber  are  cycles  per  thousand 
km  and  cycles  per  km,  respectively.  The  observed  spectral  peak  in  each  frequency 
band  is  marked  by  a  star. 
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Figure  A-3a:  The  frequency  spectrum  of  eastward-going  motions  (solid  line),  westward¬ 
going  motions  (dash-dotted  line)  and  stationary  waves  (dotted  line).  The  annual 
peak  is  most  conspicuous  in  the  spectrum  of  stationary  motions.  Figure  A-3b  The 
wavenumber  spectrum  of  eastward-going  motions  (solid  line)  and  westward-going  mo¬ 
tions  (dash-dotted  line). 
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Figure  A-4:  Frequency-zonal  wave  number  spectra.  The  peak  in  each  frequency  band 
in  Fig.  A-3  is  indicated  by  a  dot  here.  The  magnitude  of  each  peak  is  denoted  by 
the  size  of  the  circle  around  that  peak.  The  dash-dot  and  solid  curves  denote  the 
barotropic  and  first  baroclinic  Rossby  wave  dispersion  curves  with  l  =  0,  respec¬ 
tively.  At  the  lowest  frequencies,  the  barotropic  curve  is  visually  indistinguishable 
from  k  =  0.  Uncertainty  estimates  are  derived  from  the  95%  confidence  interval 
(Capon  and  Goodman,  1970)  for  the  peak  values.  The  zonal  wavenumber  interval 
represents  those  wavenumbers  where  the  background  energy  surrounding  the  apparent 
peak  is  larger  than  or  equal  to  the  lower  limit  of  the  confidence  interval. 
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Figure  A-6:  Temporally  averaged  energy  magnitude. 
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Figure  A-8:  The  temporal  evolution  of  the  energy  magnitude  of  motions  with  the 
characteristic  zonal  scale  730  km. 
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